Untitled
Untitled
Untitled
Boundary Element
Research
Edited by C. A. Brebbia
Volume 1:
Basic Principles
and Applications
With 144 Figures and 11 Tables
Springer-Verlag
Berlin Heidelberg GmbH 1984
Editor:
Dr. Carlos A. Brebbia
ISBN 978-0-387-13097-2
This work is subject to copyright. All rights are reserved, whether the whole or part ofthe material is
concerned, specifically those of translation, reprinting, re-use of illustrations, broadcasting, repro-
duction by photocopying machine or similar means, and storage in data banks. Under § 54 of the
German Copyright Law where copies are made for other than private use, a fee is payable. to
"Verwertungsgesellschaft Wort", Munich.
© Springer-Verlag Berlin Heidelberg 1984
Originally published by Springer-Verlag Berlin Heidelberg New York Tokyo in 1984
Softcover reprint of the hardcover 1st edition 1984
The use of registered names, trademarks, etc. in this publication does not imply, even in the absence
of a specific statement, that such names are exempt from the relevant protective laws and regulations
and therefore free for general use.
206113020-543210
Contributors
Carlos A. Brebbia
Editor
Preface
As the Boundary Element Method develops into a tool of engineering analysis more
effort is dedicated to studying new applications and solving different problems. This
book contains chapters on the basic principles of the technique, time dependent
problems, fluid mechanics, hydraulics, geomechanics and plate bending.
The number of non-linear and time dependent problems which have become
amenable to solution using boundary elements have induced many researchers to
investigate in depth the basis of the method. Chapter 0 of this book presents an ap-
proach based on weighted residual and error approximations, which permits easy
construction of the governing boundary integral equations. Chapter I reviews the
theoretical aspects of integral equation formulations with emphasis in their
mathematical aspects.
The analysis of time dependent problems is presented in Chap. 2 which describes
the time and space dependent integral formulation of heat conduction problems
and then proposes a numerical procedure and time marching algorithm.
Chapter 3 reviews the application of boundary elements for fracture mechanics
analysis in the presence of thermal stresses. The chapter presents numerical results
and the considerations on numerical accuracy are of interest to analysts as well as
practising engineers.
Boundary elements are becoming an accepted tool for the solution of many fluid
mechanic problems as demonstrated in Chaps. 4 and 5. Chapter 4 describes the use
of the method for solving problems in aerodynamics and hydrodynamics, free body
problems and water waves. The chapter starts with some historical remarks linking
the boundary elements to be more classical techniques such as panel methods.
Chapter 5 is dedicated to the application of boundary elements to find water waves
forces on fixed, free-floating or moored offshore structures. The authors point out the
convenience of using the simple fundamental solution of Laplace's equation to-
gether with radiation type boundary conditions, in preference to more sophisticated
fundamental solutions. The main advantage of this treatment is that it gives accurate
solutions which can be used for general types of structures without presenting any
convergence problems.
One of the advantages of boundary elements which is now beginning to be fully
exploited in commercial programs 1 is the possibility of using discontinuous ele-
ments (Chap. 6). These elements have the advantage of being able to represent dis-
continuities such as singularities better and their meshes can be refined in a much
simpler way.
The next two chapters deal with problems in geomechanics. Chapter 7 presents
the basic formulation for no-tension and problems with joints and discontinuities
such as those frequently appearing in geomechanics. The chapter also uses the
viscoplastic formulation as presented by Telles and Brebbia. Chapter 8 is ad-
dressed more specifically to the applications of the method to mining, a field in
which boundary elements have been very successfully applied. The authors present
the displacement discontinuity formulation and discuss the combination of the
method with other techniques. The chapter contains a number of interesting exam-
ples which lend weight to the authors' conclusions regarding the suitability of
boundary elements for problems in mining engineering.
Very little has been done on the topic of geometrically non-linear problems using
boundary elements. The authors of Chap. 9 present work carried out on finite de-
flections of plates using different integral formulations corresponding to diverse
plate theories. Numerical illustrations are presented for the approximate Berger
equations and their validity is discussed in this chapter. A section on extension of
the method to deal with non-linear problems in sandwhich plates and shallow shells
is also included.
The final chapter (Chap. 10) is dedicated to the interpretation of boundary in-
tegral methods using Treffiz original ideas. The chapter concentrates on the basic
mathematical definitions but is written in a way that the engineer can easily in-
terpret. The last chapter reviews the theory and in particular its applications in po-
tential and elastostatic problems, including the indirect formulations.
The book is a testimony to the strength an,d vitality of the new method which has
made considerable advances in a period of 7 years or so, giving as its 'official' date
of birth the 1st Conference in 1978 for which the name boundary elements was
used 2 and the first book published on the topic also in 1978 3 • Judging by a rule of
thumb that a method takes 15 years or so of gestation and a similar time for devel-
opment, the editor is certain that the need for this series will continue. New
volumes will continue to provide up to date information on the state of the art of the
technique and its applications at the same high standard.
Carlos A. Brebbia
../'O=~(k~)
ox ox
+ ~(k~)
oy oy
(0.3)
or an integro operator
../'0 = Sv (t-r) 0 dr (0.4)
o
which can also be written in terms of convolution notation as
L() = v (t)* () (0.5)
Most of the operators encountered in engineering formulations are differential
but for some problem areas, such as visco-elasticity and transient response, one
needs to consider integro or integro-differential operators.
Another important definition of functional analysis is the generalized inner
product. Given two vector functions u and v, the product is denoted as
<u,v)=<u·v) (0.6)
where the dot represents the scalar product of u, v at a point and the brackets may
indicate different types of operations. The most common operation is an integra-
tion over a domain Q, i.e.,
<u,v)=Ju·vdQ (0.7)
Q
and integration over a certain region of the domain. In particular, if the region is
the boundary, T, enclosing the domain, the operation is denoted by
<u,v)r=Ju·vdT (0.9)
r
Having defined the operator ...1'( ) as in (001), one can now introduce the concept
of its adjoint. The adjoint of an operator is similar, conceptually, to the transpose of
a matrix. Let u and v be two arbitrary functions and multiply (0.1) by v. We wish
to determine how the operations can be transferred to the v function, i.e.,
where...l'* ( ) is referred to as the adjoint of L. Forming the inner product over the
domain,
<...I'(u), v) = J...I'(u) v dQ (0.11 )
Q
and integrating by parts until all derivatives (we are talking about differential
operators) are passed to v gives
where P (u, v), a result of the partial integration, is called the "bilinear concomi-
tant". We will shortly examine the character of this term. The inner product of a
single function say <f, 1) is written as <I) for convenience, when no confusion as
to meaning arises.
Example 1. Consider the following second order operator;
d 2u du
...I'(u) = ao dx 2+ a] dx + a2 u . (a)
X2( d u 2 du
<...I'(u), v) = J ao -d 2 + a] - + a2 u ) v dx (b)
X, x dx
and integrating, one finds
+[v {ao~+a]u}lX2
dx x,
(c)
[t follows that
d2 d
J'* (u) = dx 2 (ao u) - dx (a] u) + a2 u
(e)
du d
P(u,v)=aov dx -u dx (aov)+a]uv
(a)
X2 x2l d 2u d 2v du dv ]
SvJ'(u)dx=S bO-d 2-d2- b ]-d -d +b2uv dx
x. x. x x x x
l{
dx
2
o - U)
+ V - d (bd
dx 2
du } +dv- ( -b o -d 2U)]
+b]-
dx dx dx 2
(b)
This result corresponds to the point where we have converted the integral such that
all the terms have, individually, the same order derivatives in x for u and v.
Continuing for two more integrations gives the final result,
l{
+ V - d (bd
dx
2
o - u)
dx 2
du } +dv- ( -b o -d 2u)]
+b]-
dx dx dx 2
-
l{
U - d (bd
dx
2
o - V)
dx 2
dV}
+b]-
dx
+du- ( -b o -d 2V)]
dx dx2
(c)
where f* () =J'(), i.e. the operator is scI/adjoint, and the R (u, v) term has the
form
2 du } +dv- ( -b o -d 2U)
R(u,v)=v {- d (bd o - u) +b]- (e)
dx dx 2 dx dx dx 2
It is conventional to define the forms of F and G after the first phase of the inte-
gration. Notice that the order of F and G is equal to half the order of J( ). Hence
is J() is of order 2 n, F and G will have terms of order n. F will involve
derivatives up to order n -I according to our convention for defining F (i.e. result
after the first n integrations). The boundary conditions are defined as
Essential
F (u) = f (a prescri bed vector)
(b)
Natural
G(u) = g (a prescribed vector).
and
When the operator is self-adjoint, one can transform the differential equation and
boundary conditions to a variational statement. We just outline the procedure here
to stress its relationship to integral equation methods, but our primary interest is in
weighted residual techniques which allow for non-self-adjoint operators, including
non-linear cases.
Consider the following differential equation,
...t'(u) + b = 0 inQ (0.17)
and boundary conditions,
F(u) = f onr l
(0.18)
G(u) = g onr2
We mUltiply equation (0.17) by v and integrate by parts until we have balanced
derivatives of u and v. This yields
S[...t'(u)+b]vdQ=SF(v)·G(u)dr+S[bv+.@(u)·@(v)]dQ (0.19)
where @() is a vector containing derivatives of u or v. We can also write (0.19) in
terms of scalar product notation,
<...t'(u) + b, v) = <F(v), G(u)r + <b, v) + <.@(u), @(v) (0.20)
As an illustration, the equation discussed in examples 2 and 3 has the following
@ ( ) operator,
2
@() = {b1l2 d 0 . b l12 dO bl12()} (0.21)
o dx2 ' I I dx' 2
I
H(u)=-b o -
2 dx2
(d2u)2
--bl -
2 dx
I (d U )2 +-b
I
2
2u
2
(0.25)
on the boundary. Since F and G are linear in u, the variational operator can be
shifted inside the bracket, i.e., 0 [F (u)] = F (ou), and the result is
o(F(u), G(u»r= (F(ou), G(u»r+ (F(u), G(ou»r (0.27)
where now F(ou) is not constrained (i.e. F(ou) =1= 0 on T,). Defining the func-
tional 1 as
12 = (H(u) + b u) + «F(u) - 1), G(U»rl + (F (u), g)r 2 (0.35)
it follows that the statement 0/2 = 0 for arbitrary ou is equivalent to the complete
set of conditions,
feU) + b = 0 in Q
F(u)=f onT, (0.36)
G(u) = g on T2
This approach to incorporating all the boundary conditions is followed in
establishing the general weighted residual form. The important difference is that
BOUNDARY INTEGRAL FORMULATIONS 7
while the above approach is only valid for self-adjoint operators, the weighted
residual technique allows us to produce boundary element type statements for any
problem.
Example 4. Consider the operator of Example 2.
F(U)={U,~~}
where
Let us now consider the set of equations represented by (0.36). Our strategy is to
propose an approximate solution a containing unknown parameters and to
establish appropriate values for these parameters. To avoid proliferation of nota-
tion, we shall drop the superscript tilde on a, and one must keep in mind that u is
now only an approximation. The errors for the respective equations are
G[ = f (u) + b =1= 0 in Q
III = F(u) - f=l= 0 on TI (0.37)
1l2=G(U)- g=l=O onT2
These errors are reduced by distributing them with prescribed weighting functions
in such a way that the product of the residuals by the weighting functions is set to
zero over Q and T, i.e.
S G[ w dQ + S III . WI dT + S 112 . W2 dT = 0 (0.38)
Q F, F2
8 BOUNDARY IN1EGRAL FORMULATIONS
where w, WI and W2 are weighting functions (the last two are vectors). By applying
(0.38) one can determine the unknown parameters in the approximate function
chosen for u.
The key issue is the choice of weighting functions and here our approach can be
guided by the result obtained with the variational method. Referring back to
section 3, we established a functional whose stationary requirement was equivalent
to the satisfaction of the system given by equation (0.36). Examination of equation
(0.34) shows that the right hand side is similar in form to (0.38) and suggests that
we take the weighting function for the self-adjoint case to be
W = ()u
WI = G(w)
w2=-F(w)
With this choice, the weighted residual expression is identical to the stationary
requirement. Then, for the self-adjoint case, we require,
for this choice, the interior integral involving u vanishes and one has only to work
with an expansion for u on the boundary. Interior node points are not required. If,
as usual in boundary elements, we choose a function such that
(0.44)
BOUNDARY INlEGRAL FORMULATIONS 9
where {) is a delta function applied at a point P and r is the distance from the
origin, the domain integral becomes the value of the function u at the point P, i.e.
Comparing (0.47) with (0.50) shows that both forms are identical when we select
W, and W2 as follows:
WI = G*(w)
(0.51)
W2= - F*(w)
or equivalently
Jew dQ + J81 . G* (w) dr- J82· F* (w) dr= 0 (0.53)
!J r, r.
Equation (0.53) reduces to (0.40) for the self-adjoint case.
10 BOUNDARY IN1EGRAL FORMULATIONS
We consider first the case where the governing equation is Poisson's equation. The
weighted residuals, are
Q
S(b w+ u'1 2 w) dQ + S (w q - u
G
~w) dr+ GS (w q - u ~w) dr= 0
~ ~
(0.57)
This latter form is the one we employ in the boundary element method, along with
the fundamental solution
(0.58)
where
I I
w=--= for 3-D
4nr 4nlx-xpl
(0.59)
I 1
w = - -In r = - -In I x - x I for 2-D
2n 2n p
Substituting (0.58) into (0.57) gives for an interior point P the well known result,
- up + S b w dQ + S (-ou w - U-
ow) dr
Q r, on on
+ S (q
r2
w- u ~w)
n
dr = 0 (0.60)
Next we consider the equilibrium equations for an elastic solid which can be
written as
I: = diva + b = .u + u + (A. + .u) V div u + b (0.61)
.u, A. are the
Lame constants; div is the divergence operator; u is the displacement
vector; a the stress tensor; V the gradient operator; b is the body force vector; and I:
BOUNDARY INTEGRAL FORMULATIONS 11
We can apply this statement in boundary elements together with the fundamental
(Kelvin) solution,
,uLlw+ (A + ,u)V divw+;; = 0 (0.65)
where;; indicates a vector of dirac delta functions, each component corresponding
to a particular direction. Hence equation (0.64) becomes for an internal point p
- up + Jb . w dQ + J (p . w - U • q) dr + J (p . w - U • q) dr = 0 (0.66)
F F. F2
This chapter demonstrated how one can establish the correct boundary conditions
for a given problem using the concepts of functional analysis.
In addition, the importance and convenience of weighted residual techniques for
formulating the governing integral statements are discussed. In particular these
techniques are useful when confronted with new problems for which an integral
statement is not available. Such is the case with many time dependent and non-
linear problems. Weighted residuals provide a simple and very convenient way to
initiate the boundary element method.
Some examples are presented in the chapter to illustrate how the governing
equations for the boundary element technique can be obtained and the meaning of
the different terms.
Bibliography
Brebbia, C. A, The Boundary Element Method for Engineers. Pentech Press, London,
Halstead Press, NY, 1978, Second Edition 1980
Brebbia, C. A, S. Walker, Boundary Element Techniques in Engineering. Butterworths,
London, Boston, 1980
12 BOUNDARY INTEGRAL FORMULATIONS
Brebbia, C. A., 1. Telles, L. Wrobel, Boundary Element Techniques - Theory and Applications
in Engineering. Springer-Verlag, Berlin, NY, 1984
Brebbia, C. A. (Ed.), Progress in Boundary Element Methods, Vol. 1. Pentech Press, London
and Halstead Press, NY, 1981
Brebbia, C. A (Ed.), Progress in Boundary Element Methods, Vol. 2. Pentech Press, London
and Springer-Verlag, NY, 1983
Brebbia, C. A (Ed.), Recent Advances in Boundary Elements. Proceedings of the 1st Int. Conf.
on BEM Pentech Press, London, 1978
Brebbia, C. A (Ed.), New Developments in Boundary Element Methods. Proceedings of the
2nd Int. Conf. on BEM, Southampton 1980. CML Publications, Southampton 1980.
Second Edition 1983
Brebbia, C. A (Ed.), Boundary Element Methods. Proceedings of the 3rd Int. Conf. on BEM,
California 1981, Springer-Verlag, Berlin, NY, 1981
Brebbia, C. A (Ed.), Boundary Element Methods in Engineering. Proceedings of the 4th Int.
Conf. on BEM, Southampton 1982. Springer-Verlag, Berlin, NY, 1982
Brebbia, C. A (Ed.), Boundary Element Methods. Proceedings of the 5th Int. Conf. on BEM,
Hiroshima, 1983, Springer-Verlag, Berlin, NY, 1983
Chapter 1
by M A. J aswon
Green's formula on the boundary [19]. Much of the mathematical foundation for
Rizzo's formulation had already been laid down by Kupradze [20]. However his
indirect approach, i.e. employing hypothetical source densities on the boundary,
did not make the same impact as Rizzo's direct approach which involved only the
quantities of engineering interest, i.e., the boundary tractions and displacements.
The essential equivalence between the direct and indirect formulations was de-
monstrated by Jaswon and Symm [21]. Rizzo's equations were solved numerically
by Rizzo et al. [22, 23] and vigorously applied by Cruse [24, 25, 26] to problems of
technological interest.
The engineering viewpoint has been fostered particularly by Brebbia in his
books [27, 28] and reports of conference proceedings [29, 30, 31] under the title
Boundary Element Method (BEM). This terminology of course brings out the fact
that useful common ground exists between BEM and the longer established Finite
Element Method (FEM). Ambitious attempts at combining BEM with FEM have
been proposed by Zienkiewicz [32] and by Brebbia and Georgiou [33]. Although
these developments are encouraging, it must be said that a practically effective
numerical and error analysis of the discretisation procedure over curved bound-
aries remains to be completed. Some important theoretical contributions by
Wendland [34] point the way forward.
The basic scalar and vector formulations are reviewed in this chapter together
with the necessary background of theory. Some previously unpublished material is
included.
Superposing the contributions from all over DB, we construct the simple-layer
potential
V(p) = S g (p, q) a(q) dq; pc B i , Be, DB (1.2)
oB
where B; denotes the interior domain enclosed by DB and Be denotes the infinite
domain exterior to DB. Although the main properties of V are intuitively clear on
physical grounds, its precise mathematical properties depend upon the smoothness
of DB and of a as will now be discussed.
A REVIEW OF 1HE THEORY 15
4) V(p) exists at every pc oB, and it is continuous at p with respect to its neigh-
bouring values in B i , Be i.e.
showing that
og·
J-'dq=O (1.18)
BB or
A REVIEW OF 11IE THEORY 17
aVe o r r age
-;)- = -;- Jg e (p, q) (Jo dq = (Jo J : ; - dq, etc. (1.21)
~ ~M u~
aVe)
(-;)-
ur r=a
= - 4n (Jo =V;(p) , (7OV..)
ur r=a
=0= VEep) (1.22)
I
4) W(p) exists at every pc oB, but it jumps on passing from oB to neighbouring
points outside oB. More precisely, if Pi, Pe are points on each side of oB near p,
then lim W(Pi) = W(p) + 2 n P (p)
Pi~p ; pc oB (1.28)
lim W (Pe) = W (p) - 2 n P (p)
Pe ...... P
I
oW(p) og
-ot- = o~ at
(p, q)i p (q) dq; pc oB (1.30)
As will be seen later, these simple topological results play an important role in the
theory of Fredholm integral equations. For ease of reference, we note the cor-
responding results:
Jg(p,q)~dq=-4n; pcBi (1.36)
aB
= - 2n; pc oB (1.37)
0; pc Be (1.38)
These may be verified for the sphere of radius a by choosing
p= (0,0, z), q = (r sin () cos '1', r sin () sin '1', r cos ())r=a,
so that
dq = a 2 sin () d() d'l' ,
g (p, q) = (z2 + r2 - 2z r cos ())~~/;,
Green's Formula
It has been noted above that the simple-layer and double-layer potentials are
harmonic functions under broad conditions. However, an arbitrary harmonic func-
tion may not always be representable by such potentials. For instance, the
harmonic function rp = k (a cons.) can not be represented by a simpre-Iayer
potential inside the unit circle. Again, the harmonic function rp = r- 1 can not be
represented by a double-layer potential in the infinite domain exterior to a closed
surface. To construct a more powerful potential representation for harmonic func-
tions, we posit a harmonic function rp in B i , which assumes a continuous set of
boundary values rp (q), and boundary normal derivatives rpi (q), as q runs over oB.
These boundary data, regarded as source densities, generate the double-layer
potential
J
g (p, q)i rp (q) dq , (1.39)
aB
and the simple-layer potential
- J g (p, q) rp~ (q) dq , (l.40)
aB
which have the properties associated with W, V respectively for HOlder-continuous
boundary data on a Liapunov surface oB. Superposing (1.39) and (1.40) yields the
identity
Jg (p, q)i rp(q) dq - Jg (p, q) rpj(q) dq = 4n rp(p); PCBi (1.41)
oB oB
valid for any harmonic (jJ in B i . This is Green's formula [35].
Formula (1.41) can be readily extended from pc Bi to pc oB. Thus (l.40)
remains continuous, but (1.39) jumps by - 2 n rp (p), at oB so providing the
boundary formula
J g(p,q)jrp(q)dq- J g(p,q)rpj(q)dq=2nrp(p); pcoB (1.42)
oB oB
20 A REVIEW OF 1HE 1HEORY
This differs essentially from (1.41) in that rp on the right-hand side is now a
boundary value of rp, i.e. from the same set as enters into the first integral, i.e.
(1.42) is a functional relation between rp, rpj on vB which ensures their compatibil-
ity as boundary data. As p crosses from vB to Be, the double-layer integral suffers
a second jump - 2 n rp(p) yielding Green's reciprocal theorem
S g (p, q); rp(q) dq - S g (p, q) rpi(q) dq = 0; pc Be (1.43)
aB aB
as follows directly from the Gauss divergence theorem. Starting from (1.43), and
reversing our steps, we successively recover (1.42), (1.41). Of course, this procedure
hinges upon the validity of the boundary jumps (1.28), which can only be justified
by the same kind of limiting analysis as would be involved in the direct proof of
(1.41).
Green's formula may be readily adapted to a regular harmonic function I in
Be. which assumes boundary values I(q), and boundary normal derivatives I: (q),
as qruns over vB. Corresponding with (1.41), (1.42), (1.43) we have
S g(p,q)~/(q)dq- J g(p, q) I: (q) dq=4nl(p); peBe (1.45)
aB aB
=2n/(p); pevB (1.46)
=0; pcBi (1.47)
respectively, where
1
I=--Ipl-l J/:(q)dq+Olpl-2 as Ipl-+oo (1.48)
4n aB
All these exterior formulae are important in their own right, and (1.47) may be
exploited to justify the single-potential representations (1.2), (1.24). Thus super-
posing (1.41) and (1.47), assuming rp given and I arbitrary, we obtain the
generalised representation
J g (p, q)i[rp(q) - I(q)] dq - J g (p, q) [rpj (q) + I: (q)] dq = 4n rp(p); pcB;. (1.49)
aB aB
Two natural possibilities for I now arise:
1) 1= rp on vB, giving the representation
- Jg (p, q) [rpi( q) + I; (q)] dq = 4 n rp (p) ; pcB; , (1.50)
aB
Simple-Layer Formulations
capacitance problem of electrostatics. In this case the charge density A satisfies the
equation
Jg (p, q) A(q) dq = I; P e oB (1.57)
oB
and we note that
I) A > 0 on oB, so providing the capacitance
x= S A(q) dq > 0; (1.58)
oB
2) A generates the simple-layer potential
S g (p, q) A (q) dq = 1, pcB;; (1.59)
oB
3) A satisfies, in addition to (1.57), the normal derivative equation (1.65) presented
below.
Equation (1.57) has been solved numerically for various closed surfaces as a
means of computing their electrostatic capacitance.
According to the interior Neumann existence theorem, there exists a unique (up
to an arbitrary constant) harmonic function rp in B;, which assumes prescribed
continuous normal-derivative values rp; on oB subject to the Gauss condition
S rpi(q) dq = 0 (1.60)
oB
To construct rp in B;, we introduce the representation (1.54) as before, with rp;
assumed Holder-continuous on a Liapunov surface oB. This yields the boundary
relation
J
g;(p,q)a(q)dq-2na(p)=rp;(p); pc oB (1.61)
oB
which may be viewed as a Fredholm integral equation of the second kind for a in
terms of rpi. Despite the singularity in the kernel, see (1.16), classical Fredholm
theory applies: the homogeneous equation
S g;(p,q)a(q)dq-2na(p)=0; peoB (1.62)
oB
has a non-trivial solution a = A, corresponding with the non-trivial solution f1 = 1 of
the adjoint homogeneous equation
Jg (p, q); f1 (q) dq - 2 n f1 (p) = 0; p coB (1.63)
oB
as may be seen from (1.34), so that a solution exists subject to the orthogonality
condition
J
oB
J
f1 (p) rpi(p) dp = ifJ;(p) dp = 0
oB
(1.64)
i.e. the Gauss condition (1.60). By suitable normalisation, A may be identified with
the unique solution, A, of equation (1.57), since this yields the normal derivative
equation
J
g;(p,q)A(q)dq-2nA(p)=0; peoB (1.65)
oB
Accordingly, equation (1.61) has the general solution
a= ao + k A (1.66)
A REVIEW OF mE THEORY 23
The same procedure applied to (1.61) gives the expected result S rp; (p) dp = O. An
important companion result is aB
S qJ (pH (p) dp = S a (p) dp (1. 74)
aB aB
as may be proved by operating with S ... A (p) dp upon both sides if (1.55) and
~gain invoking Fubini's theorem: aB
Double-Layer Formulations
The representation
qJ (p) = S g (p, q); ,u (q) dq ; pcB; (1. 75)
aB
provides a classically preferred alternative to (1.54) for Dirichlet problems, where
,u appears as a hypothetical Holder-continuous source density to be determined. In
1
principle ,u = - (rp - f), but of course f is not available ab initio. An effective
4n
way forward is to note that qJ remains continuous at iJB, whereas the integral jumps
by - 2 n,u (p) at p c iJB, so yielding the boundary relation
S Je(p) {qJ(P) __
ICI} dp= S Je(p)qJ(p)dp-c S JeI(PI) dp=O (1.84)
iJB P iJB BB P
This value of c ensures that equation (1.83) has a general solution fl = flo + k,
where flo is a particular solution and k is an arbitrary constant, which generates the
unique potential
J g(P,q);flo(q)dq+k J g(p,q);dq; pcB e (1.85)
BB BB
I.e.
J g(p, q); flo (q) dq; pc Be (1.86)
iJB
Direct Formulations
Both the Dirichlet and Neumann problems may be formulated directly through
Green's boundary formula (1.42). Thus, given 91 on oB (interior Dirichlet
problem), it becomes a Fredholm integral equation of the first kind for 91; in terms
of 91, i.e.
Jg(p,q)qJ;(q)dq= J g(p,q);qJ(q)dq-2nqJ(p); pcoB (1.87)
BB iJB
This has a unique solution as already noted in connection with equation (1.55).
With 91, 91; both known on oB, we may generate 91 throughout Bi via Green's
formula (1.41). It requires proof, however, that 91; satisfies the Gauss condition
(1.60): operating with J... Je(p) dp upon both sides of (1.87), and interchanging
iJB
the order of integration where appropriate, we obtain
l.e.,
S rp;(cOdq S g(q,p)A.(p)dp= Srp(cOdq S g;(q,p)A.(p)dp-27O S rp(p)A.(p)dp
oR . oR oR oR oR
i.e.,
S rp; (cO dq = 270 S rp (cO l (cO dq - 270 S rp (p)A. (p) dp = 0 (1.88)
oR oR oR
on bearing in mind (1.57), (1.65) with p, q interchanged.
Given rp; on oB (interior Neumann problem), equation (1.87) now reads as
S g(p, cO;rp(cO dq-27Orp(p) = S g(p,cOrp;(cOdq; pcoB (1.89)
oR oR
i.e. a Fredholm integral equation of the second kind for rp in terms of rp;. This has
an associated homogeneous equation
S g(p, cO; rp(cO dq - 2n rp(p) = 0; pc iJB (1.90)
oR
and adjoint equation
S g;(p,cOl(cOdq-27Ol(p)=0; pcoB (1.91 )
oR
which are mathematically equivalent to (1.80), (1.81). Consequently, equation
(1.89) only has a solution subject to the orthogonality condition
where g' (Pa q~) signifies the traction component in the IX-direction at p generated
by a unit point-force acting in the n-direction at q. Clearly column I defines the
traction vector at p generated by a unit point-force acting in the I-direction at q,
etc. Finally, corresponding with g(p, q)' we construct the traction dyadic
where row I defines the traction vector at q generated by a unit point-force acting
in the I-direction at p, etc. It may be readily proved that column I defines a
singular displacement vector at p, i.e., that generated by a unit traction-source
associated with the I-direction at q, etc., in line with the fact that g(p, q)' may
function as a unit dipole-potential generated at q. By the same token, row I of
(1.103) defines the singular displacement vector at q generated by a unit traction-
source associated with the I-direction at p. Any individual component of (1.103) or
(1.104) carries two possible interpretations: either the traction component gen-
erated by a unit point-force or the displacement component generated by a unit
traction-source. The interpretation will always be clear from the context. We note
that g' (p, q) stands for g~ (p, q) or gj (p, q) as the case may be, and similarly for
g(p, q)'.
The simple-source density a now becomes a vector simple-source density
a = <a], a2, a3), so allowing us to regard (1.2) as a vector simple-layer potential
with components
Va(P) = J g(Paq~)a~(q)dq; 1X,I7=I,2,3 (1.105)
aB
This has properties at oB entirely analogous to those of the scalar simple-source
potential, e.g. formulae (1.12), (1.13) may be read as traction formulae, and it
defines an eIastostatic displacement field for any choice of p. These properties
have been proved by Kupradze for the linear isotropic elastic continuum, but we
may conjecture that they also hold for the general linear anisotropic elastic
continuum. Similarly, the double-source density /1 becomes a vector double-source
density /1 = </1]' /12, /13), so allowing us to regard (1.24) as a vector double-layer
potential with components
(1.106)
A REVIEW OF lHE lHEORY 29
Somig/iana's Formula
Green's formula (1.41) now reads as Somigliana's formula [37], i.e. it represents an
arbitrary displacement vector qJ as the superposition of a vector simple-layer
potential and a vector double-layer potential, generated respectively by the
boundary tractions and boundary displacements associated with qJ. Green's
boundary formula (1.42) now reads as Somigliana's boundary formula, which
provides a vector functional relation between tractions and displacements on oB.
Green's reciprocal theorem (1.43) now reads at Betti's reciprocal theorem.
Corresponding exterior formulae hold for a displacement field which remains
regular at infinity, so allowing us to introduce the generalised Somigliana formula
(1.49). The boundary-conditions of linear elastostatics are all covered by (1.97),
reading qJ, qJ', y as vectors and rI., fJ as scalars. Assuming the validity of the
fundamental existence-uniqueness theorems of linear elastostatics, we may readily
prove the validity of the vector representations V, W with a, J.1 identified by the
vector equations (1.51), (1.53) respectively.
the point-force at P, and that (2) their resultant moment about any axis balances
its moment about this axis. Easy deductions from (1.109) are
S g(p, q)i fls (q) dq = 2n fls (p); s=I,2, ... ,6, pcvB (1.110)
oB
=0; pc Be (1.111)
in parallel with (1.34), (1.35) respectively.
Kupradze Formulations
Kupradze was the first to propose that the vector interior Neumann problem could
be formulated by the vector integral equation (1.61), with the associated homo-
geneous systems (1.62), (1.63). Naively applying classical Fredholm theory to this
system, it follows that a solution only exists if
S !p~(p) fls(P) dp = 0; s = 1,2, ... ,6 (1.116)
oB
i.e. if the prescribed tractions form a self-equilibrated distribution of forces and
moments over vB, as expected on physical grounds. Subject to (1.116), the general
A REVIEW OF THE THEORY 31
where (10 is a particular solution and ks are arbitrary scalar coefficients. This
generates the class of displacement fields
(1.122)
Jg (p, q)~,u (q) dq + 2 n,u (p) = IP (p) - a . g (q, p)q=o - b /\ V . g (q, p)q=o (1.129)
as
always has a solution, of the form (1.122), which generates IP in Be via (1.125).
Rizzo Formulations
As first proposed by Rizzo, Somigliana's boundary formula (1.42) may be
exploited to attack vector Dirichlet and Neumann problems. Thus, given the
displacement vector qJ on oB (interior Dirichlet problem), equation (1.87) now
reads as a vector integral equation of the first kind for lPi in terms of qJ, This has a
unique solution lPi, as already noted for the vector equation (1.55). Operating with
J... As(p) dp upon both sides of (1.87), we see that lPi satisfies the relations
as
J lPi(p),us(p)dp=O; s=I,2, ... ,6 (1.130)
as
as expected on physical grounds,
Given the traction vector lPi on oB (interior Neumann problem), equation (1.89)
now reads as a vector integral equation of the second kind for qJ in terms of lPi.
Subject to condition (1.130), this has a solution which appears as
6
These satisfy the relations (1.7) - (1.15) adapted as appropriate, except that (1.9)
gets replaced by
V(p) = log Ipi J a(q) dq -lpl-2 J (p. q) a(q) dq + 0 Ipl-2 as Ipi -> 00 (1.135)
08 08
These satisfy the relations (1.25)-(1.31) adapted as appropriate, except that (1.27)
gets replaced by
W(p)=0Ipl-2 as Ipl->oo (1.138)
A convenient feature of two-dimensional theory is the Cauchy-Riemann relation
log Ip - q Ii = logi Iq - pi = - aq
iJ()
(q - p) ; pcB; (1.139)
34 A REVIEW OF lHE lHEORY
rp~
instance: given. = K (a cons.) on r = a i.e. (~~) = K, there exists a unique
harmomc functIOn r~a
rp = K a log r; r ~ a (1.146)
conforming to (1.143). By contrast, however, no solution exists subject to bounded
behaviour at infinity. The difficulty is surmounted in classical theory by imposing
the exterior Gauss condition J
rp~ (q) dq = 0. This forces the choice K = 0, which is
88
consistent with the class of solutions
rp(r) = c (any constant including zero); r ~ a (1.147)
One of these solutions, i.e. rp = 0, is uniquely determined in the class (1.146) by
putting K = 0.
A REVIEW OF lHE lHEORY 35
T-Contours
A distinguishing feature of two-dimensional theory is the existence of contours,
termed T-contours, for which the equation
J 10glp-qlA.(q)dq=l; peoB (1.148)
oB
does not have a solution. To identify these contours, start with a contour oB
defined by the Cartesian equation f(x, y) = 0, and construct the family of similar
contours f(xlm, ylm) = 0 where m > 0 is a continuously varying parameter. Then,
according to theorems proved by Jaswon & Symm:
(i) one, and only one, member ofthis family is aT-contour;
(ii) on T, and only on T, the homogeneous equation
J log I p - q I A. (q) dq = 0 (1.149)
oB
has a non-trivial solution.
To locate T, we label the m th contour oBm and introduce the normal derivative
equation
J
log; p - q I A. (q) dq + n A. (p) = 0; p e 0B (1.150)
oB.
on oBI (= oB). This has a non-trivial solution A corresponding with the non-trivial
solutionfl = I of the adjoint equation
J10glp-qlifl(q)dq+nfl(p)=0; peoB, (1.151)
oB.
Clearly A generates the potential
J 10glp-qIA(q)dq=k(acons.); peoB (1.152)
As already indicated in connection with (1.144), the unit circle is a T-contour for
the family of concentric circles: by symmetry equation (1.148) has the solution
36 A REVIEW OF lHE lHEORY
or (ii) J qJ~ (p) dp '* 0, which implies J 0" (q) dq =1= 0 and therefore the represen-
oB oB
tation (1.156).
Plane Elastostatics
Vector potential theory may be adapted to two dimensions by using appropriate
expansions for the dyadics g (p, q) or g' (p, q) or g (p, q)'. A simpler alternative
approach is to work with a scalar biharmonic function X (V2 V2 X = 0), which
signifies either the Airy stress function (plane stress or plane strain) or the
transverse deflection of a thin plate. According to a theorem of Hadamard, X is
uniquely determined in B; if x, xi are given on iJB - data readily available from
the boundary tractions or boundary deflections as the case may be. To construct X
in B;, we employ the Almansi representation
x=r2 tp+lf/; r2=x 2+y2, v 2tp=O, V21f/=O (1.169)
with tp, If/ represented as simple-layer potentials. Thus, writing
where (UI, U2, 0) signifies the displacement vector and tP is the harmonic conjugate
of IP. The formulae (1.173) may be proved directly, but they also follow from the
Papkovich-Neuber representation
Given UI, U2 on vB, formulae (1.173) provide a pair of coupled integral equations
for IP (and therefore also tP) and 'II on vB, so enabling IP, tP, 'II and therefore also
UI, U2 to be generated in B i . Details are given in a recent paper by Symm and
Bhattacharyya [38] building upon earlier work by Jaswon and Symm. Very recently
[39] the representation (1.172) has been employed to attack traction problems as an
alternative to (1.168).
Uniqueness for exterior problems requires X = 0 (r) as r -> 00. There is no
difficulty in adapting the representation (1.169) to this requirement: we write
X = r21P + 'II + c (1.175)
where IP. 'II are provided by (1.170) and c is a constant balanced by S a (q) dq = 0
aB
(which ensures IP=O(r- l ) and therefore x=O(r) as r-> 00). However, the
adaption of (1.172) to exterior domains presents considerable difficulties which go
beyond the scope of this review.
References
I Fredholm, I., Sur une c1asse d'equations fonctionelles. Acta Math. 27,365 - 390 (1903)
2 Hess, J. L., Smith, A M. 0., Calculation of non-lifting potential flow about arbitrary
three-dimensional bodies. Report No. E. S. 40622, Douglas Aircraft Co., Long Beach 1962
3 Hess, J. L., Smith, A M. 0., Calculation of potential flow about arbitrary bodies. In:
Progress in Aeronautical Sciences, Vol. 8. (D. Kuchemann, Ed.). Pergamon Press, London
1967
4 Jaswon, M. A, Ponter, A R. S., An integral equation solution of the torsion problem.
Proc. Roy. Soc. A, 273, 237-246 (1963)
5 Jaswon, M. A, Integral equation methods in potential theory, I, Proc. Roy. Soc., A, 275,
23- 32 (1963)
6 Volterra, v., Theory of Functionals and of Integral and Integro-Differential Equations.
Dover, New York 1959
* These formulae refer to plane strain, but they also hold for plane stress on replacing v by
v/(l + v)
A REVIEW OF THE THEORY 39
7 Symm, G. T., Integral equation methods in potential theory, II. Proc. Roy. Soc., A, 275,
33-46 (1963)
8a Symm, G. T., External thermal resistance of buried cables and troughs. Proc. LE.E.E.,
116 (10),1695-1698 (1969)
8b Symm, G. T., Computation of potential in a multiwise proportional counter of arbitrary
cross-section. Nucl. Instrum. Methods, 118,605-607 (1974)
9 Symm, G. T., Capacitance of coaxial lines with steps and tapers. In: Recent Advances in
B.E.M, Proc. pI Int. Conj. B.E.M, Southampton Univ., (c. A Brebbia, Ed.). Pentech
Press, London 1978
10 Symm, G. T., An integral equation method in conformal mapping. Num. Math., 9,
250- 258 (1966)
11 Symm, G. T., Numerical mapping of exterior domains. Num. Math., 10,437-445 (1967)
12 Symm, G. T., Conformal mapping of doubly-connected domains. Num. Math., 13,
448-457 (1969)
13 Papamichael, N., Whiteman, 1. R., A numerical conformal transformation method for
harmonic mixed boundary value problems in polynomial domains. 1. App. Math. Phys.
(ZAMP), 24,304- 316 (1973)
14 Gaier, D., Integralgleichungen erster Art und konforme Abbildung. Math. Zeit., 147,
113-129 (1976)
15 Hough, D. M., Papamichael, N., An integral equation method for the numerical
conformal mapping of interior, exterior and doubly connected domains. Num. Math. (in
the press)
16 Jaswon, M. A, Maitit, M., Symm, G. T., Numerical biharmonic analysis and some
applications. Int. 1. Solids Structures, 3,309-332 (1967)
17 Jaswon, M. A, Maiti, M., An integral equation formulation of plate bending problems.
1. Eng. Math., 2 (I), 83-93 (1968)
18 Rizzo, F. 1., An integral equation approach to boundary value problems of classical
elastostatics. Quart App. Math., 25 (I), 83 - 95 (1967)
19 Jaswon, M. A, Some theoretical aspects of boundary integral equations. App!. Math.
Modell., 5,409-413 (1981)
20 Kupradze, V. D., Potential Methods in the Theory of Elasticity. Israel Program for
Scientific Translations. Jerusalem 1965
21 Jaswon, M. A, Symm, G. T., Integral Equation Methods in Potential Theory and
Elastostatics. Academic Press, London 1977
22 Rizzo, F. 1., Shippy, D. 1., A Method for Stress Determination in Plane Anisotropic Elastic
Bodies. 1. Composite Materials, 4, 36-60 (1970)
23 Vogel, S. M., Rizzo, F. 1., An integral equation formulation of three dimensional
anisotropic elastostatic boundary value problems. J. Elast. 3 (3), 203 - 216 (1973)
24 Cruse, T. A, Numerical solutions in three dimensional elastostatics. Int. 1. of Solids
Structures 5, 1259-1274 (1969)
25 Cruse, T. A, An application of the boundary-integral equation method to three-
dimensional stress analysis. Computers & Structures 3, 509-527 (1973)
26 Cruse, T. A, An improved boundary-integral equation method for three-dimensional
elastic stress analysis. Computers & Structures 4, 741-754 (1974)
27 Brebbia, C. A, The Boundary Element Methodfor Engineers. Pentech Press, London 1978
28 Brebbia, C. A, Walker, S., Boundary Element Techniques in Engineering. Butterworth
1979
29 Brebbia, C. A (Ed.), Recent Advances in B.E.M. Proc. ]SI Int. Con! B.E.M. Southampton
Univ., Pentech Press, London 1978
30 Brebbia, C. A (Ed.), New Developments in B.E.M., Proc. 2nd Int. Can! B.E.M. C.M.L
Publications, Southampton 1980
31 Brebbia, S. A (Ed.), B.E.M., Proc. 3rd Int. Seminar. Irving, California. C.M.L. Publica-
tions, Springer 1981
32 Zinkiewicz, O. K, Marriage a la Mode-Finite Element and Boundary Integral Method.
In: International Symposium on Innovative Numerical Analysis in Applied Engineering
Science, Versailles. CETIM Publications, Paris 1977
40 A REVIEW OF lliE 1HEORY
2.1 Introduction
~ = k V 2u in Q xT (2.1)
at
satisfying the initial condition
u(x, to) = uo(x) in Q v (2.2)
and appropriate boundary conditions. In heat conduction problems we assume the
following forms for the boundary conditions:
- Dirichlet boundary condition:
u(x,t)=u(x,t) on rjxT (2.3 a)
- Neumann boundary condition:
q(x, t) = q(x, t) on r2 x T (2.3 b)
42 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
is a particular solution of (2.4). Notice that x E Rd and all t > t' and r represents
the euclidean distance between points x and x', i.e.,
r=(r·r)1/2, r=x-x' (2.6)
It can be shown (see reference [I]) that this function has the following properties:
lim u* (x, t; x', t') = b(x - x') (2.7 a)
t'--+t
where 150 stands for the Dirac distribution centered at the origin;
S u* (x, t; x', t') dQ (x) = S u* (x, t; x', t') dQ (x') = I (2.7b)
Rd Rd
The function u (x, t; x'; t') is called the fundamental solution of equation (2.4). It
represents the field of temperature produced by an instantaneous point source of
heat placed at point x' and instant t'. Next we demonstrate how this fundamental
solution can be used to obtain an integral representation for the solution of the
heat equation. Consider the following
ou
u=u(x,t'), fj(=kV 2u, (2.8a)
ou*
u* = u* (x, t; x', t'), --= - kV2 U * (2.8b)
ot' '
from which we see that u and u* are solutions of two adjoint equations. We also
known that:
o ou ou*
- ( u u*) = u* - + u - - = k (u*V2 u - uV2 u*) (2.9)
ot' ot' ot'
APPLICATIONS IN lRANSIENT HEAT CONDUCTION 43
Integrating both sides of this equation over the cylinder r x (to, t) we obtain
t o t
I I - , (uu*)dQ(x')dt'=kI I(u*V 2U-uV 2u*)dQ(x')dt' (2.10)
toQ ot to Q
The left hand side of this expression can be modified by interchanging the order of
integration, i.e.,
/ 0 / 0
J I;;; (u u*) dQ (x') dt' = I I ;;; (u u*) dt' dQ (x') = I [u u*lio dQ (x')
/0 Q vt Q to vt Q
= I u (x, t) u* (x, t; x', t) dQ (x') - Ju (x, to) u* (x, t; x', to) dQ (x')
Q Q
If the point x lies in the interior of Q then, by the well known properties of the
Dirac distribution, c(x) = 1. If x is on the boundary r of Q then c(x) equals the
fraction of solid angle subtended by r at x, relative to the solid angle of the sphere
in Rd [21]. For smooth points on r, c(x) = 112.
Applying Green's second identity to the right hand side of (2.10) we find that
ou ou* )
I (u* V 2u - u V2 u*) dQ (x') = I ( u* - - u - - dr (x') (2.13)
Q r ~ ~
Collecting this result and (2.12) and abbreviating slightly the notation, we arrive at
the expression t
cu=Suou~dQ+S S(u*q-uq*)dt'dr (2.14)
Q r to
This relation shows that the value of the function u at interior points (c = I) and
for any instants t > to can be explicitly evaluated by integration, once the initial
temperature field and boundary values of temperature and flux are known. Or, to
use other words, (2.14) gives an integral representation for the solutions of the heat
equation.
We notice that in (2.14) the fundamental solution flux is
ou* 1
q*=k--=- r·n·u* (2.15)
on 2 (t - t')
We have two kinds of space integrals in the representation (2.14), one in the
domain Q and another on the boundary r. Now we can show that, if Uo is a
harmonic function in Q, the space integral can be transformed into a boundary
integral. Indeed, applying Green's second identity to Uo and another function rp we
have (orp ou )
1(uoV2rp-rpV2Uo)dQ=J uoa;;-rp ono dr (2.16)
44 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
But, by definition,
(2.17)
J
AUo uo dQ = (u o !: - ~:o) dr
cp (2.19)
The integration of (2.18) for two and three dimensional problems brings no special
difficulties. For example, for the two dimensional case employing polar co-
ordinates we have
J...~ (r
r dr
drp)
dr
I
4nk(t-to)
exp [
- r2
4k(t-to)
1 (2.20 a)
1 1
rp= - 4 n E] (4kt- to») - 4 n In (4k(;~ to») - :n
1 (2.22)
cp = - 4 n e] (4k(;- to) )
equation, i.e.,
-
ou = k V2 u in T
Q x (2.24a)
ot
subject to the following two types of boundary conditions,
u= a on rlxT (2.24b)
q= q on r2xT (2.24c)
where a is equal to the prescribed temperature values on r l and q represents the
prescribed flux at points on the boundary r 2. As before r I and r 2 are
complementary parts of the boundary r of the domain Q. We recall expression
(2.14) and, in view of (2.24 b, c), rewrite it as
1
We note that the operator :Yfotransforms functions over r 2 into functions over r
and the operator ~ transforms functions over r I into functions over r. Then, the
expression (2.25) can be simply written as
(2.27)
Equation (2.25), or its equivalent (2.27), is also valid at points on the boundary.
Making x E r l we obtain the following relation
:Yfou2-~ql=b-Cii (2.28 a)
and making x E r 2 we similarly find the following expression
C U2 +:Yfou2 -~ql = b (2.28b)
Defining a new operator ;;y;> and a new function b such that,
;;Y;>U2 = C2 U2 + :Yfou2 (2.29 a)
and
b= b- CI a (2.29 b)
with CI and C2 functions on r given by
CI = C, Cz = 0 for x E r l
(2.30)
c]=O, C2=C for XEr 2
46 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
--tn
where the symbols L,I and L,2 stand for summation over the nodes of fl and of r2.
The parameters U2i~ and qlia represent the nodal values of U2 and ql at node ion r
and node IY. on T. Since by definition the lfIi and ({J~ are interpolation functions, at
the nodes Xi and ta we have that
If/i(Xj) = {)ij and qJa (tp) = {)aP (2.33)
Introducing relations (2.32) into equation (2.31) we obtain
L,2 L, (:Yr' lfIi ({Ja) U2ia + L, I L,;:9' (Ifli ({Ja) qlia = b (2.34)
To determine the nodal values U2ia and qlia we collocate equations (2.34) at points
Xj on r and points !p in the interval T. Defining
which, once solved, yields the nodal values U2ia and qlia'
---i
I ~ I
I I I
i I I I
I I I I
_+-!__---+!_ _----!I_ _-,1-!_ _---'--____
fo f1 12 fn_1 fn ·f
Fig. 2.2. Piecewise constant elements in time
where
[HII] = Hjli!. [G)l] = Gjlil and {bj} = bjl (2.40 a)
[HII] is a N(r) x N(r 2) matrix, [GIl] is a N(r) x N(r l ) matrix and {b} a N(r)
vector. Since N(r) = N(r l ) + N(r 2) we have the same number of equations as
unknowns. We remember that at each point of r either the temperature or the flux
(but not both) is an unknown.
The integrals involved in matrix [HII] are
I,
[HHl == Hjli! ==(:YI'!lfi'fJa)(Xj, tl) = S S q*(Xj' tl; x', t') !lfi(X/) dt' dT(X/) (2.40 b)
t2 10
[G)l] == Gjlil == (~ !lfi 'fJ!J.)(xj. tl) = S SU* (Xj, tl; x', t') !lfi (x') dt' dr (x') (2.41)
t,lo
To obtain {b l } we need to compute boundary integrals that have the same form
as (2.40, 41) plus an integral over the domain.
Recalling expressions (2.5, 15) we can write that
I, 1
[HJl]=- S S , r·n!lfiu*dt'dr(x / ) (2.42 a)
t, 10 2(tl-t)
and grouping all the terms depending on the variable t' we obtain
The integration on time can be done analytically. In fact making the change of
variable
?
s (2.43)
4k(t1 - t') ,
= 4k
,2
exp(-so)=~exp(-~)
,2 4kLlt
(2.44)
[GJ}] ==
1,
S S 4n
1
k(t _ t') exp -
l 4k(t _ t')
,2 1'!Iidt'dr(x')
t,10 i l
(2.47)
where CI) 1
El (s) = S- exp(- z) dz (2.4S)
s z
is the exponential integral. Approximations suitable for computer applications can
be found in reference [S, ch.5]. Introducing (2.47) into (2.46) we have that
1
[GJ}]=- S E\ (-r2-) '!Iidr(X') (2.49)
4nk t, 4kLlt
To evaluate these integrals we proceed as before, except that now the integrand
develops a singularity of the type In, as , --+ o. This means that numerical
integration over the elements containing the collocation point has to take this into
consideration. According to (2.21) we can write
E 1 (s) = - In s + el (s) (2.50)
where el (s) is a smooth function. Then the integral in expression (2.49) can be split
into a regular part and a singular part with a logarithmic singularity. We can apply
the above considerations for the regular integral. For the singular one must employ
special integration rules such as those due to Berthod-Zaborowski [9]. The whole
50 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
done once at the first time-step, saved and used in the forward and backward
substitutions required at all later time-steps.
In this procedure we require the computations of domain integrals JU~ u~ dQ
(j
even if the inital temperature field Uo is zero. To evaluate those integrals two
methods can be enployed. One, which we designate by Method A, is to compute U a
at internal points coinciding with the integration points of the cells. The other,
Method B, is to compute U a at the nodes of the mesh and construct a finite element
type approximation. As the number of integration points usually is greater than the
number of nodes this second variant requires evaluation of the temperature of a
smaller number of internal points, thus being more rapid and demanding less
computer memory. The first method however tends to be more precise, but if the
size of the cells is too small some integration points can become too close to the
collocation points where the fundamental solutions changes rapidly and this may
offset the greater theoretical accuracy.
Procedure ll. This procedure consists in starting all time integrations from the
initial instant to up to the current instant tn. Suppose the unknown nodal vectors
{u~} and {q1} have been calculated at t a , rI. = I, ... , n - I. Equation (2.36) allows
one to advance the solution to time tn. Considering that the collocation time tp in
(2.36) equals tn, we obtain,
n
L, ([Hna] {un - [Gna] {t]1}) = {b n} (2.52)
a~l
Isolating the unknowns {uq} and {tl]} on the left hand-side we have the following
system of linear equations
n-l
[Hnn] {U1} - [G nn ] {q'l} = {b n} - L, ([Hna] {un + [Gna] {tlm (2.53)
a~l
which, once solved, delivers {U2} and {tln. This procedure requires the formation
of a triangular array of matrices that can be depicted as
/3=1
/3=2
/3=3
/3=n
and similarly for [GPaj. We recall that tp is the collocation time and the ta are the
nodes in the interval T = (to, t n ). It can be easily shown that the matrices connected
by the diagonal lines in this table are equal. Therefore at step n only the matrices
[Hnlj and [Gn1j and the right-hand side vector {b n} have to be calculated anew.
Albeit this procedure is more elaborate than Procedure I it presents some
advantages. One is that computation of values of temperature at internal points is
not required. The other rests on the fact that to evaluate the domain integral
52 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
JUo uti dQ in (2.26c) we need only to discretize the support of the function u at
Q
time to. When Q is an unbounded domain this can be a major benefit. Also if the
initial temperature field Uo (x) is harmonic in Q the above integral can be
transformed into a boundary one, as discussed in section 2.2. It should be remarked
that the left-hand side of (2.53) is identical with that obtained with Procedure I.
Therefore the LU factorization needs to be performed only once and saved. At
each time step the solution of (2.53) requires only· a forward and a backward
substitution phases.
There is so far no general proof that the above system is solvable. However the
special cases treated by [6] and [14] and the numerical experience gained from
BEM practice validate the approach.
As we have remarked before, during the time marching process one system of
the type (2.54) has to be solved for each time step using the same matrices but
with different right-hand sides. Thus the LU factorization method seems the most
adequate. As soon as the matrices [H! - G] are formed they are factorized and this
factorization saved. At later steps only the forward and backward substitution
phases need be performed. The matrix [H! - G] is a dense matrix, hence the
number of operations (multiplications/divisions) that the factorization involves is
N (t)3 /3, the forward and backward substitution phases requiring N (t)212 each
[15, ch.3].
We have been using subroutine SGECO of Linpack [16] to estimate the
condition number of the resulting matrices for several test cases found that in all of
them the condition number was typically less than 102, indicating good stability for
the system (2. 54 b).
~
::::>
a a 1.0
:,;
0.
0.4
E
~
5.0
4 6 10
b c t-
Fig. 2.3. a Geometry. b Discretization. c Results for r = 1
By a simple change of variable this problem can be transformed into one with
zero initial condition, thus making Procedure II of section 3 specially suitable.
The mesh employed is depicted in Figure 2.3b. We see that because of
symmetry one can discretize only one quarter of the hole circumference. Only 2
quadratic elements and 5 nodes were required to obtain the results shown above in
Figure 2.3 c.
The numerical integration were performed with 4 Gauss-Legendre points for the
regular integrals and with 3 special points [9] for those with logarithmic
singularities. The time step employed was L1l = 0.5.
The results obtained for the temperature at r = I are shown in Figure 2.3 c and
they agree well with those of references [6] and [17].
Example 2. This is one-dimensional test problem used in reference [18]. The region
Q is the interval 0 < x < I subject to the initial condition Uo (x) = sin n x and to
the Dirichlet boundary conditions u (0, t) = u (1, t) = 0 for all t > O. The analytical
solution is known to be u (x, t) = exp (- n2 t) sin n x.
The mesh employed for this example is shown in Figure 2.4 Procedure I,
Method A and Procedure I, Method B yield similar accuracy, Method A proving to
be marginally better.
With linear elements the best accuracy was obtained with a time step
L1 t = om and with quadratic elements with L1 l = 0.005. The accuracy of
Procedure I, Method A remains almost unchanged for a large range of time steps.
Procedure I, Method B is more sensitive to the choice of time step. Typical results
are presented in Figure 2.4.
Example 3. The aim of this example is to show the applicability of the BEM to
more complicated cases. The problem is to find the temperature field in a gas
turbine rotor blade with cooling holes. The convection coefficient h was obtained
[com [19, 20]. The geometry and mesh employed are depicted in Figure 2.5 as well
as some isothermals for several times. A total of 45 quadratic elements and 90
nodes were used. In order to draw the isotherms the temperature was also
computed at 30 interior points. Integrations were performed as in Example I.
54 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
1.0,-----------------, 1.0..----------------,
- analytical - - analytical
• t= 0.01 o BEM-A
o t = 0.05
0.8 0.8
~ 0.6
:::J
EOJ
Cl.
~0.4
Example 4. This example consits on a section of a structure with a hole filled with
air. The geometry is shown in Figure 2.6. The boundary conditions are of mixed
type, the top and bottom surfaces are radiating according to the following formula:
. (t - 4)
top surface: q= 900 sm n -I -6 -W/m 2 (t in hours); absorptivity = 0.7
. (t - 5)
bottom surface: q= 83 sm n -I-6 - W 1m2 (t in hours); absorptivity = 0.5
On the vertical surfaces the radiation intensities vary linearly on time, as follows.
4 OW/m2 0 OW/m 2
lO 400 W/m2 10 100 W/m 2
ILl 200W/m 2 15 600W/m 2
20 OW/m 2 20 OW/m 2
APPLICATIONS IN TRANSIENT HEAT CONDUCTION 55
Holes T(°el H! W
/cm1k)
1 480 C.067
2 500 Q.062
3 470 1.069
c 4 460 J074
r
"'"
~E
3:
en
"
.<;:
40 35
~ 30
23/
25
~
21 r,\
\. 21,
c ~o 25
23,
~o/
45 ___
40
35
~
25
23_
30
~
21-
fV
J'
2V
2h
d ~o 25 ~5
40
35
30 25
23
21
On the lateral surfaces the absorptivity is 0.5. The thermal properties of the section
are: conductivity = 1.4 W1m· K, density = 2400 kg/m 3, specific heat = 960 Jlkg . K.
The ambient air temperature also has a linear variation on time,
t(hours)
4
15
20
The initial condition, assumed at t = 8 hours, is a uniform temperature u = 20°C.
The mesh consists of 70 linear elements and 78 boundary nodes. The number of
internal points required to draw the isotherms was 90. The time step was chosen
L1 t = I hour. The temperature of the air in the hole was computed via a energy
balance equation, that is, the energy lost by convection through the sides of the
hole was used to raise the internal energy of the air. The results obtained using
Procedure II in this complex example are shown in Figure 2.6.
Conclusions
The examples discussed in this chapter and which range from simple tests to more
complex practical applications demonstrate that the BEM can be easily applicable
to the solution of linear heat conduction problems. The so-called Procedure II
allows for the discretization only on the boundary thus leading to an effective
reduction on the dimensionality of the problem. In the cases where comparisons
where possible the BEM yielded good accuracy with time steps which are
substantially larger than those employed in the finite element or finite difference
methods. This conclusion has been confirmed by other authors, specially Wrobel
and Brebbia in reference [6]
It is important to point out that the application of the BEM to the heat equation
still requires more research work, namely with respect to the three dimensional
case. Also it should be stressed that more effort is needed to develop computa-
tionally efficient codes.
References
1 H. S. Carslaw,1. C. Jaeger, Conduction of Heat in Solids. Oxford U. Press 1959
2 B. A Boley, 1. H. Weiner, Theory of Thermal Stresses. 1. Wiley 1960
3 F. 1. Rizzo, D. H. Shippy, A method of solution for certain problems of transient heat
conduction. AIAA 1.,8,11,2004- 2009 (1970)
4 Y. P. Chang, C. S. Kang, D. 1. Chen, The use of fundamental Green's functions for the
solution of problems of heat conduction in anisotropic media. lnt. 1. Heat Mass Transfer,
16,.1905-1918 (1973)
5 R. P. Shaw, An integral approach to diffusion. lnt. 1. Heat Mass Transfer, 17,693-699
(1974)
6 L. C. Wrobel, C. A Brebbia, Time dependent potential problems. Chapter 6 in Progress
in Boundary Element Methods. Vol. 1. (Ed. C. A Brebbia), Pentech Press, London,
Halstead Press, USA 1981
58 APPLICATIONS IN TRANSIENT HEAT CONDUCTION
Notation
In this chapter the following general notation has been adopted:
the d-dimensional euclidean space.
an open connected set of Rd representing the domain occupied by the body.
the boundary of Q, which is supposed to be sufficiently smooth in the sense
of Kellog [21].
r ,r2
J complementary parts of r.
T= (to, tf) an interval of time, were to represents the initial instant and tf the final instant.
x, x' coordinates of points in Q or on r.
t, t' instants in T.
u (x, t) a real function over Q x T representing the temperature of the body.
uo(x) a real function over Q representing the (initial) temperature at instant to.
o
ot (-) partial derivative with respect to time t.
3.1 Introduction
Thermal stress analysis is one of the most important subjects in engineering and
technology. Several attempts have been so far made for analysis of thermoelastic
problems in steady heat conduction states by the boundary element method (Rizzo
and Shippy, 1978; Danson, 1981, for example). However, as far as the authors are
aware, there are few investigations for thermal stress analysis in transient heat con-
duction states, although the transient heat conduction problems are often studied
by means of the boundary element method (Wrobel and Brebbia, 1981). This is
partly because the boundary integral formulation includes inevitably domain in-
tegrals which arise due to the temperature distributions in space and Jime, and in
the transient heat conduction states they can not be transformed into the
corresponding boundary integrals. Therefore, we have to innovate a more efficient
numerical implementation for such problems.
In this paper, a relatively simple boundary element formulation and its numeri-
cal implementation are presented for the two-dimensional thermoelastic problems.
The proposed method of solution makes full use of the numerical techniques de-
veloped in the finite element methods for evaluating the values of stresses and
strains at internal points and also domain integrals.
A couple of sample problems in transient heat conduction as well as thermo-
elasticity are computed to demonstrate the potentiality of the proposed method.
Then we apply the method to the stress intensity factor calculation of cracked
structural components subjected to different thermal shocks. Finally, we discuss
briefly the numerical accuracy of the proposed method.
Boundary r = r; + I1
Fig.3.1. Coordinate system and notation
tive. It is assumed throughout in this study that the material constants are invari-
ant in space and time. The boundary conditions for the heat conduction equation
(3.1) are such that
() = e on r 1
q == o(J/on = q on r2 (3.2)
q = h «(Ja - (J) on r3
where r = r 1 + r 2 + r 3, (Ja denotes the ambient temperature and h the heat trans-
fer coefficient. The superimposed bar denotes a prescribed value and 0 (')/on the
directional derivative along the outward unit normal. In addition, the initial con-
dition at the time t = 0 is expressible as
(3.3)
We now introduce the fundamental solution (J* to the heat conduction equation
(3.1). This solution should satisfy in the infinite domain of the same properties as
those of the problem under consideration the following differential equation:
x 'i]2 (J* + (J* + £5 (x, x') £5 (t, r) = 0 (3.4)
in which £5{-) denotes the Dirac delta function. We also introduce the heat flux q*
which can be defined as
q* = o()* / on. (3.5)
Now, we consider a weighted residual statement of the governing equation (3.1)
in which the fundamental solution is used as the weight function. It can be given
by t
S J(x'i]2()-O) ()* dQ dr= O. (3.6)
OD
From the above expression, we can derive in the usual manner (Brebbia and
Walker, 1980) the following boundary integral equation:
t t
C () + x JJq* () dr dr = x SS ()* q dr dr + S ()* (Jo dQ (3.7)
or or D
FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES 61
in which C is the coefficient that depends only on the geometrical properties of the
boundary. For a smooth boundary point C = 0.5, and for an internal point C = 1.
The fundamental solution e* for a two dimensional isotropic medium can be
expressed as (Carslaw and Jaeger, 1978)
t-
and
It should be noted that the time-dependent fundamental solution has the following
property:
()* (t, r) == 0 for t ~ r (3.10)
which represents the principle of causality. If the boundary integral equation (3.7)
is solved for the unknowns () and q on the boundary, we can obtain the tem-
perature change at an arbitrary point throughout in the domain by using the in-
tegral equation (3.7) with the coefficient C = 1.
Now, we can calculate thermal stresses in the elastic domain under considera-
tion, having known the temperature distribution over the whole domain of interest.
As has been well known, the governing differential equations of the thermoelastic
problems can be expressed as
Dijkl Uk,lj (x) + eij (),j (x) = 0 (3.11)
where D ijkl and eij are the well known material constants of the linear elastic body
(Malvern, 1969). For the so called Hookean elastic solids D ijkl are expressible in
terms of the shear modulus G and Poisson's ratio v, and the components of e ij are
expressed in terms of G and the coefficient of linear thermal expansion IY..
The boundary conditions for Eq. (3.11) can be prescribed in such a way that
(3.12)
where (Jij denotes the stress tensor and Ui the displacement vector.
We now introduce Kelvin's fundamental solution Uti to the governing equation
(3.11), which should satisfy the differential equation:
D lamb Utl,ab (x, x') + 8km b (x, x') = 0 (3.13)
in the infinite domain of the same linear elastic properties as in the problem to be
analyzed. Here we denote by 8km the Kronecker delta.
We also consider a weighted residual statement of the governing equation (3.11)
as in the heat conduction problems. That is,
S(D lamb Um,ab + ela (),a) Uti dQ = O. (3.14)
Q
Making use of the Gaussian divergence theorem and taking into account the prop-
erties of the Kelvin's fundamental solution, we can obtain the integral equation
which holds at an arbitrary internal point x in the domain Q. Then, taking the
62FRAC1URE MECHANICS APPLICATION IN THERMOELASTIC STATES
limiting process in which the point x approaches the boundary point y and evalu-
ating the singular integrals, we can arrive at the following boundary integral
equation:
and
U:I =---l--l
8 n G (1- v)
(3 - 4 v) In (~) bkl + _o_r__o_r ]
r OXk OXI
(3.17)
=- 1
4 n (1- v) r
lonor {
- -
(1-2v)bkl+2----
OXk OXI
or or}
where P:I denote the boundary tractions corresponding to the fundamental solu-
tion U:I.
Ii
j-l j
The time integration in Eq. (3.20) can be analytically carried out, and this yields
S q* (t;, r)
I.
dr= -- 1- ( -or ) [exp (- ak-I) - exp (- ak)] (3.21 b)
1.-1 2 n r on %
where r2
ak=---- (3.22)
4% (t; - tk)
The column vector {w} in the above equation represents the contributions of the
initial condition (3.3). The coefficient matrices [H] and [G) have a triangular form
due to the principle of causality, and can be expressed as
1
::.O::L
o 0 ... 0
(3.24)
~gfl
gll
o
g22
0 ... 0
0 ... 0
~
[G) = (3.25)
For the solution of Eq. (3.23), we can proceed successively from the first time
step. Namely, we first solve the system of equations:
(3.26)
for the boundary unknowns of {8 (II)} and {q(ll)}' The column vector {w (II)} de-
notes the contribution of the initial condition at the first time step. Once the
boundary unknowns at the first time step have been determined, we can proceed to
the second time step. For an i-th time step, we have
- ----+-I-
[h u] {8 (Ii)} = [gii] {q (Ii)} + {tV (Ii)} + L ([gik] {q (lk)} - [h ik ] {8 (lk)}). (3.27)
k=1
Since {8(lk)} and {q(tk)} for k=1,2, ... , i - l are known from the solutions at
previous time steps, we have only to obtain the boundary unknowns at the i-th
time step.
It should be noted that a six-point Gaussian numerical quadrature can be used
for computing the components of the matrices [H) and [G) except the diagonal
ones. For evaluating the diagonal components of the sub-matrices [h ii ] and [gii]'
special care must be taken. The diagonal terms of [h ii ] has the singularity of 1Ir,
but they vanish due to the orthogonality of rand n, which makes or/an = O. The
diagonal terms of [gii] have a logarithmic singularity which is directly integrable.
As shown in Wrobel and Brebbia (1981), expanding the exponential-integral func-
tion in series, we can evaluate the diagonal components [gii]jj in the following
closed form:
[gii]jj=JL{2- Y-In(aj )
2n
+ I
m= I
(_l)m-I a? ,}
m (2 m + 1) m.
(3.28)
where
17
a. = ] (3.29)
J 4X(li- li_l)
in which y is Euler's constant, i.e. y= 0.57721566 ... and Ij is the length of thej-th
boundary element.
FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES 65
After having known the temperature distribution in space and time by means of
the boundary element analysis mentioned above, we can calculate the domain-
integral term OJk in Eq. (3.15). This is done by subdividing the internal domain Q
into small cells and using the standard Gaussian numerical quadrature rule. Then,
the usual boundary element method available for elastostatic problems can be
applied to the solution of the thermoelastic problems under consideration. Using
the same discretization scheme as in the heat conduction problems, we can obtain
the following system oflinear equations:
where [A] and [D] are coefficient matrices which can be calculated from the fun-
damental solutions (3.17) and (3.18), respectively. In equation (3.30) the column
vectors {u} and {p} denote the nodal displacement and nodal traction vectors,
respectively, and {j} the column vector corresponding to the temperature distri-
bution throughout in the domain.
Equation (3.30) is now solved for the unknown components of {u} and {p} by
taking into account the boundary conditions (3.12). Then, we can compute the
thermoelastic stresses and strains produced at an arbitrary time. Applying this
solution procedure to every time step, we can determine the whole behavior of the
thermoelastic problems.
When we are concerned with fracture mechanics application of the above solu-
tion procedure, the method of sub-regions is very useful, in particular for the stress
intensity factor computation of a crack with arbitrary shape and location.
We now consider a two dimensional domain divided into two sub-regions as
shown in Fig. 3.3. The subdivision is made in such a way that the interface rI
includes the cracked surface. The outer boundaries of the sub-regions (Domain 1
and Domain 2) are denoted by r] and r2, respectively. Denoting by subscript I the
variables on the interface, we can express the discretized set of Eq. (3.30) for the
two sub-regions. Namely, we have for the sub-region I
(3.31)
(3.32)
where {ail and {pi} are the nodal displacement and traction vectors on the bound-
ary T i , respectively.
The nodal -displacements {an and nodal tractions {pH on the interface T[ must
satisfy both the compatibility and equilibrium conditions. These conditions can be
expressed as
{aJ} = {an == {all
(3.33)
{p}} = - {pi} == {h}
Making use of Eq. (3.33) for Eqs. (3.31) and (3.32) and rearranging them, we can
obtain
(3.34)
On the cracked surface, however, the compatibility condition (3.33) can not be
used. Since in most cases the cracked surface is traction-free, instead of Eq. (3.33)
the following conditions are applicable to the cracked surface:
Substitution of Eq. (3.35) into Eq. (3.34) leads to the global set of equations for the
unknowns on the outer boundary Tl + T2 and also on r i . Instead of a direct appli-
cation of the Gaussian elimination to the solution of these equations, the following
elimination procedure can be recommended for practical use. Namely, we first
eliminate the supplemented unknowns on the interface T I , and then the real
unknowns on the outer boundary Tl + T 2 . If this stepwise elimination procedure is
used in the method of sub-regions, the system matrix to be inverted has the same
order as that of the original boundary element equations without sub-regions.
When only smaller computers are available, the so-called block elimination proce-
dure, which is one of the standard numerical techniques in the finite element
methods, can also be applied to the solution of the boundary element equations
associated with the method of sub-regions. It is, however, beyond the scope of this
article to describe the details of this procedure.
FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES 67
'""f~"'C
0 BEM
x, - Exaet
100
'C 'C
80 80
<V
::::>
-0
~ 60
E
w
t-
40
20 20
o 4 8 em 10 0 8 em 10
x- x-
Fig. 3.4. Time variations of temperature Fig. 3.5. Time vanatlOns of temperature
for case of one-side cooling for case of two-side heating
68 FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES
Crack tip Xl
Fig. 3.6. Coordinate system at crack tip
We now define the coordinate system at the crack tip as shown in Fig. 3.6. Al-
though we may apply different procedures to the evaluation of the stress intensity
factors, we shall use in this study the following two methods. The first one is called
the stress method in which the values of stress intensity factors are estimated by
linear extrapolation using the stress values in the vicinity of the crack tip. When
our attention is paid to the stresses 0"22 and 0"12 on the axis Xl (If! = 0), we can ob-
tain K] and Kn as the limiting values defined by
K] =lim
r-+ 0
[~(0"22)'I'=0]
.
1 (3.36)
Kn = lim [~ (0"\2)'1'=0]
r .... O
The other one is named the displacement method in which the crack surface dis-
placements in the vicinity of crack tip (If! = n) are used to obtain K] and Kn by
t1n
linear extrapolation, i.e.
2G
K] = --lim - (U2) 'I' = 1t
l+Jc r .... O r
(3.37)
0.6 r---~--------,
c;::- ',«-
'"b-
0.3 "
0 .....
0.2 L...--_-'---_-'--_-'--_--L.._-.-l
o 0.02 0.04 0.06 0.08 0.1
r/b-
Fig. 3.7. Circular cylinder with two opposite Fig. 3.8. Example of linear extrapolation
radial cracks
FRACTURE MECHANICS APPLICATION IN TIIERMOELASTIC STATES 69
1.0
0.8 0
0 0
0
r 0.6r- 0
0
• • • •
c;:: O.4r- • •
•
0
0
•
0.2
1.0,-----------------,
0.8
r·
b
6
o'C 0.4
..;: 0 BEM
-FEM
-100'C 0.2
Fig. 3.10. Edge-cracked plate subjected to Fig. 3.11. Normalized stress intensity factor
linear temperature gradient FJ = K I Ko of edge-cracked plate
alb
BEM FEM
0.8.-------------,
0.6
20 40 s 60
Time
t (sec) KIIKo
0.5 o.5,--------------,
0.4 0.4
r 0.3
..:::
!
ct'
0.3
0.2 0.2
• Cooling shock
o Heating shock
0.1 1 0.11~0-'1~~~~~~10~~~~10'1~~s~10J
10- 10
Time Time
Fig. 3.13. Effects of alb on Fl' R;lRo = 0.8, Fig. 3.14. Effects of boundary conditions
on inner surface () = - 100 °C, on outer on Fl' R;lRo = 0.8, alb = 0.2 and heat flux is
surface () = 0 °C and over crack surface heat zero over crack surface; x, on inner surface
flux is zero () = - 100°C and outer surface () = 0 °C;
0, on inner surface () = 0 °C and outer sur-
face () = 100°C
steady heat conduction states, we first calculate the stress intensity factors under
mechanical loadings and then proceed to the cases of thermal loadings.
Now, we show the results obtained for the rectangular plate with an oblique
edge crack as shown in Fig. 3.15.
In Fig. 3.16 comparison is made between the results obtained by the present
BEM and Freese's analytical ones (Sih, 1973) for the case of the plate subjected to
FRACTURE MECHANICS APPLICATION IN TIfERMOELASTIC STATES 73
2.0
•
z
To ('[)
1
1.8
1.6
.1·
/
1
Ll2 1.4
Y
1.2
<l:'
1.0
0.8 0.7
Ll2 ~
J
Q
0.6 "c(/ 0.6
"';,e
0-- ....
",.0:/ 0.5 r
== ;~}Freese
% .~
"
Fig. 3.15. Rectangular plate with oblique 0.4 cE'
edge crack
0.3
Fig. 3.16. Values of F, and F" for rect- • P =30' o p= 45'
angular plate with oblique edge crack sub-
jected to uniaxial tension a 0.1 0.2 OJ 0.4 0.5 0.6
a/b-
uniaxial tensile stress 0'0 at its ends. The plate aspect ratio is assumed as Lib = 2.
The values of K( and K" are normalized as
(3.38)
There are no appreciable differences between both the results. They are also sum-
marized in Table 3.5.
We now discuss the case where the plate is initially kept at uniform temperature
To (0C) and suddenly at time t = 0 its one side is cooled into O°C as shown in
Fig. 3.15. It is assmued that Lib = 4, b = 10 cm and the plate is in a plane strain
state. The stress intensity factors are norm<:.lized in this example such that
F,=K/1Ko, F"=K,,IKo, Ko =E (J.. To yrna/(l- v). (3.39)
74 FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES
1.2.-------------,
0.8 r - - - - - - - - - - - - - - - - ,
1.0
alb = 0.1
0.5 0.8
0.5 o BEM
~
Lt:' 0.4
0.2
O~~ __ ~ __ ~~~~-L---L-L~
Fig. 3.17. Time variations of F j and F" for Fig.3.1S. Values of F j and F" for
rectangular plate with oblique edge crack sub- rectangular plate with oblique edge
jected to thermal loading crack in steady thermal state
Table 3.6. Time variations of stress intensity factors in rectangular plate with oblique edge
crack
t (sec) alb
F] Fu F] Fn F] Fu
Figure 3.17 shows the time variations of F j and Fn for various aspect ratios when
qJ = 45
0 • The results are also summarized in Table 3.6. In these examples a steady
state is approached in about 500 seconds. In Fig. 3.18 comparison is made between
the present results and Nishitani's ones (Nishitani, 1975) for the steady state heat
conduction.
Finally, we compute the center-cracked rectangular plate as shown in Fig. 3.19.
It is assumed that the plate is initially kept at uniform temperature O°C and
suddenly at time t = 0 both the sides are heated into 10 (0 C). The plate is in a
plane strain state. The stress intensity factors are made nondimensional in the same
manner as in Eg. (3.39).
The results obtained are shown in Fig. 3.20 for the case qJ = 0 0 , and also sum-
marized in Table 3.7 for various values of qJ and a/b.
FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES 75
0.4r---------------,
O'C I L/2
0.3
~n;;;i;gl
To ('C)
heating
-
0.1
L/2
~b-~ O~~--L-~~-LL,~~~-L~
6 10 10 2 S 10 3
Time
Fig. 3.19. Center-cracked rect- Fig. 3.20. Time variations of FI in case of
angular plate rp = 0° for center-cracked rectangular plate
Table 3.7. Center-cracked rectangular plate subjected to heating shock from opposite two
sides
t (sec) rp= 0° rp = 45 0
FI FH FI FH FI FH
10 0.301 0.288 0.263 0.270 0.236 0.249 0.312 0.246 0.325
20 0.338 0.306 0.254 0.296 0.255 0.276 0.320 0.248 0.321
30 0.317 0.288 0.230 0.277 0.230 0.250 0.288 0.221 0.302
50 0.255 0.224 0.180 0.211 0.162 0.191 0.219 0.157 0.235
100 0.133 0.119 0.101 0.105 0.080 0.099 0.111 0.093 0.122
200 0.038 0.036 0.028 0.038 0.D25 0.036 0.036 0.028 0.038
Domain 1
0.7 , - - - - - - - - - - - - - ,
0.6
~ack --0-00..
0.576 ~,
'\
t 0.5 0.',
Domain 2
c:::'
~
,
,,
Q4 ~
Inner
boundary 0.3 '----'----'---'--'----'------'-----'
1 J.... 1 ....L 1
90 70 60 50 40
a b l/N-
Fig. 3.21. Example of boundary element dis- Fig. 3.22. Number of boundary
cretization (a) and internal cell subdivi- elements (N) versus numerical ac-
sion (b) curacy
1.1.----------------, 1.1 r - - - - - - - - - - - . . . . . ,
Exact o Exact 0
[J
1.0 -®---r---y--..L----- o- - 1.0 f-oo--a---'l....---J--------
W 6 • W
a a
OO'C x 00.C
~ 09
....... ~ 0.9
o·c{: ]100'C
:::E :::E
I..W
o 5005 =
I..W
o·c ; 100'C
=
0.8 • 3005 0.8
- Ii. 6 1005
- Ii.
0.7 '---'--_..L-_--'---....J....._ _ _- - ' - _ - - - - ' _ - ' 0.7 L...----'_---'_----'_----'_----'_----.J
1
IOif
J....
32 16
l
10"8
1 1 o 20 40 60 80 100 120
Time interval (5)
1I(Number of internal cells)
Fig. 3.23. Influence of internal cell subdivision Fig. 3.24. Influence of time stepping
on accuracy
FRACTURE MECHANICS APPLICATION IN THERMOELASTIC STATES 77
References
1. Brebbia, C. A, Walker, S., Boundary element techniques in engineering, Butterworths, Lon-
don 1980
2. Carslaw, H. S., Jaeger, J. c., Conduction of heat in solids, 2nd ed., Clarendon Press, Ox-
ford 1978
3. Danson, D. J., Boundary element methods, ed. by C. A Brebbia, pp. 105-122, Springer-
Verlag, Berlin 1981
4. Hellen, T. K, Cesari, F., Martin, A, Int. J. Pressure Vessels and Piping 10, 181-204,1982
5. Malvern, L. E., Introduction to the mechanics of a continuous medium, Prentice-Hall,
New York 1969
6. Nishitani, H., Trans. Japan Soc. Mech. Engrs., 41, 1103, 1975
7. Sih, G. C. (ed.), Mechanics of Fracture, Vol. I, p. 51, NoordhofflHolland 1973
8. Wrobel, L. c., Brebbia, C. A, Progress in boundary element methods, Vol. I, ed. by C. A
Brebbia, pp.192-212, Pentech Press, London 1981
Chapter 4
4.1 Introduction
In this chapter we make an effort to review the applications of boundary methods
to fluid mechanics. At the outset we wish to give our definition of boundary
methods. First, we exclude such techniques from the general class of finite element
methods. Although some of the language and a few of the numerical techniques of
the two methods have merged, they are historically quite separate.
However, we do classify as boundary methods a number of techniques that have
little or no difference in concept, application, or history but appear under different
names. These include: the Boundary Integral Equation Method, the Boundary
Element Method, the Boundary Integral Method, the Surface-source (or Surface-
velocity, Surface-doublet, Surface-singularity) Method, and Panel Methods *. Nor
do we wish to differentiate between the "direct method" (where the physical
variables appear as solutions of the integral equation) and the "indirect method"
(where the integral equation solution leads to the strengths of singularities from
which the physical variables can be calculated). We view all these techniques as
variations of the method of Trefftz [I].
We are hopeful that this review will show something of the wide range of
applications that boundary methods have in fluid mechanics problems, and that it
will encourage engineers not using such techniques to begin applying them. We
have attempted to cite references of historical importance and those that are recent.
Many more references can be found in the papers that are cited herein. For those
beginning to apply boundary methods, the emerging books represent the easiest
way to become acquainted with the procedures [2, 3,4, 5, 6, 7, 8].
4.2 History
Far from being a new theory, the boundary element method has a long history.
Although it became popular in the 1960s, those who developed it for use on
computers stood squarely on the shoulders of the applied mathematicians of the
* We concede that Green's function methods and panel methods, for example, are different
in their mathematical hasis and that some techniques using one method are not readily
apparent using the othu. Both, however, involve the use of fundamental singularities and
integral equations, and in ,orne cases are indistinguishable.
APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS 79
19th and early 20th centuries. The basic integral equations for potential theory
were formulated by Green in 1828 [9J.
The first application of boundary methods to fluid mechanics (and probably the
first numerical calculation using boundary method in any application) was in the
amazing paper by Trefftz in 1917 [1]. He solved the problem of a round jet issuing
from a hole in an infinite tank with the objective of calculating the contraction
coefficient. Such a calculation was not repeated until 1967 (see Hunt, [10]) using
modem computing equipment. Thus, Trefftz not only made the first boundary
integral calculation, he was first to solve an axisymmetric problem with the ()-
integration formulated in terms of elliptic integrals. He used a method of
successive approximations to satisfy the integral equation. His calculated contrac-
tion coefficient was 0.611, which was assumed correct for many years, although
Hunt [10] and others have since corrected it to 0.58.
A doubly symmetric potential flow about an elliptic cylinder was solved by
Prager in 1928 [11]. His calculation divides the boundary into elements and uses
simultaneous equations to solve the integral equation. [However, the replacement
of the integral with a finite sum and the use of simultaneous equations had also
been done by Nystrom in 1928 [12] who used the integrals of Fredholm [13] to solve
a torsion problem.] Due to the double symmetry Prager was able to obtain
reasonable results with only three algebraic equations.
Calculations using distributed singularities along the axis of two dimensional or
axisymmetric bodies became routine in the early 20th century. That method
was, however, restricted to flow around long slender objects and could unexpectedly
fail for certain shapes. To remedy this problem Lotz [14] (later Fliigge-Lotz, who
was to become a distinguished professor at Stanford University) developed the
surface-source method. She applied source distributions directly on the surface of
a body. Integration was done numerically or graphically and the integral equations
were solved by iteration. The method was later adopted by Vandrey in 1951 [15]
who used a surface distribution of sources and vortices.
Axially symmetric flow has been a frequent application of boundary methods.
Weinstein in 1948 [16] used sources which were placed on rings, disks, or cylinders
to solve this problem. Van Tuyl in 1950 [17] and Sadowsky and Sternberg in 1950
[18] continued that development. Landweber [19] used what is essentially a surface-
source method with velocity as an unknown in the integral equation. Vandrey [15],
[20] attempted to simplify the labor involved for practical calculations.
Thus, it is clear that boundary methods had become well known and were in use
during the 1950's. The labor involved in their use was an obstacle to general
popularity. The use of the digital computer, just becoming a practical tool at that
time, was a natural step which should be classified as evolutionary rather than
revolutionary. Such a use was probably first done by Smith and Pierce in 1958 [21]
although many groups, especially in the aircraft industry and at the National
Advisory Committee on Aeronautics (NACA, the predessor to NASA) had also
begun to use the computer. Having resolved the difficulty of the artihmeticallabor,
the uses of boundary methods began to expand rapidly. Three dimensional
problems in potential flow were solved by Hess and Smith in 1962 [22]. Other
notable contributions were made by Davenport in 1963 [23], Jaswon in 1963 [24],
Symm in 1963 [25], and Chaplin in 1964 [26].
80 APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS
in which Voo is the free stream velocity, rf is the strength of the ith vortex on the
body, and rr is the strength of the ith vortex in the wake. In order to have a single
valued potential, the total vorticity must be zero, i.e.,
N M
L rf+ L rr=O
i~l i~N+l
(4.4)
The rE are adjusted so as to fulfill the boundary condition alP/an = 0 on the body.
The wake vortices are those shed from the corner of the body and acquire the
strength of the corner vortices as they are shed. They are advected by a finite
difference stepping scheme
(4.5)
in which rW is the position of the vortex and u is the velocity vector which results
from the total collection of vortices and the parallel flow. Inamuro, et al. [32] found
it necessary to limit the 8-component of the core velocity of each wake vortex, i, as
follows
-1r ·W (Q ~ a)
1
UiO= 2nQ I
_Q_r'Y
(4.6)
(Q < a)
2nr?- I
82 APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS
in which Q is the distance from the vortex center and a = 2.24 Vv?f where v is the
kinematic viscosity and t1 is the time that has elapsed since the ith vortex was
shed.
Although the vortex separation method treats only inviscid flows, viscous flows
have also been calculated by boundary methods. White and Kline [32] used
boundary elements to compute flows through an axisymmetric diffuser in which
the boundary layer may be either attached or separated. Only the potential part of
the flow, however, was computed by using boundary elements. The wall geometry
was altered by adding the boundary layer displacement thickeness, b*. Separation
points were found by standard boundary layer techniques. The flow downstream of
a separation point was computed with the streamline through the separation point,
'fIs, acting as a free surface (Figure 4.2). Along 'fIs is constant, implying that the
velocity is constant. Thus, the potential is taken as
(4.7)
in which tJ>sp is the potential at the separation point, s - ssp is distance along 'fIs
from the separation point, Vsp is the velocity at the separation point (equal to the
exit velocity of the diffuser), and ssp is the position of the separation point. After an
initial guess for 'fIs, the problem is iterated whereby the normal derivatives
along 'fIs, computed from the boundary element method, are used to adjust the
position of 'fIs so that otJ>lon -4 0 along the free streamline. Details of the iteration
method are found in White and Kline [33].
In a remarkable series of papers dating from 1970, Wu and his colleagues (see
Wu [34], for a review) have applied boundary methods to viscous flow problems.
They solve the incompressible Navier-Stokes equations (Wu and Thompson [35])
ov 1
-+(v'V) v=--Vp+vV2 v (4.8)
at Q
V'v=O (4.9)
in which the symbols have their usual meaning of velocity, time, density, pressure,
and kinematic viscosity. They divide the flow into its kinematic and dynamic
components. The kinematic component is described by the continuity equation
(4.9) and the dynamic component is described by the vorticity vector
w =Vx v (4.10)
I
I
I
-
I
I
I
I
Fig. 4.2. Separation in an axisymmetric diffuser (after White and Kline, 1975)
APPLICATIONS OF BOUNDARY ELEMENT ME1HODS TO FLUID MECHANICS 83
By taking the curl of (4.8), the pressure is eliminated and thus the dynamic
component of the problem is governed by
ow = w . V v- v· V w + V V 2 w
- (4.11)
at
In the velocity and vorticity values are known at some time k, (4.11) can be used
to advance the vorticity values to time k + 1. These new values could then be used
in (4.8) and (4.9) to compute velocities at time k + 1. Instead, equations (4.8) and
(4.9) are written in integral form in which
Rainfall
Rainfall
Unconfined region
x
Fig. 4.3. A complex groundwater basin (after Liggett and Liu, 1983)
axes, and with the equations of continuity V .v = 0, (4.17) can be written simply
(Liggett and Liu, [6])
(4.18)
in which the operator V now applies to the rotated and stretched coordinate
system.
Boundary conditions for a porous media flow consist of the usual Neumann and
Dirichlet conditions except where a free surface exists. In the latter case, the
elevation, 11, of the free surface becomes a part of the problems. Then the boundary
conditions become
cfJ=11 (4.l9a)
(/J = G -a - (/J -a
S1 dto S [O(/J oG ] ds + S(/JG_I dV (4.23)
o r n n R 1=0
The foregoing equations and boundary conditions represent a "mix and match"
situation. In groundwater calculations several different such situations in which
different equations and/or boundary conditions can apply to portions of the overall
problem. These include situations which are basically two-dimensional, axisym-
metric or three-dimensional. Many such situations have been assembled into a
comprehensive program which zones in both the horizontal and vertical directions
(Liggett and Liu [6]). The zones are connected by equating pressures and flow
along the interzonal boundaries.
86 APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS
Many special situations exist in porous media calculation. For example, the
calculations in the vicinity of wells is important in most situations. For regional
calculations wells appear as singularities and may be represented as point or line
sinks. At the other extreme the detailed flow in the vicinity of a well can be treated
by the solution of an axisymmetric problem with or without a free surface
depending on the nature of the aquifer (Lennon et al. [37], [38]). The intermediate
case of a well field, in which several wells interact, is a full three-dimensional
problem (Lennon et al. [38]). In all of these cases boundary methods appear to
offer outstanding advantages over domain methods. In the first instance the
boundary method makes dealing with both the singularity and the far field
especially easy. In the last case the ease of handling the free surface and the
efficiency in three-dimensional problems become important.
One important problem that has not been solved by boundary methods is the
case of flow in unsaturated media. That problem is so highly nonlinear that it
appears boundary methods have no advantage over a domain method. In the case
of saturated-unsaturated flow, however, boundary and domain methods can be
corr bined so as to take maximum advantage of each (Liggett and Liu [6])
The complexities introduced by a free boundary are often more serious than
nonlinearities in the governing differential equations. The nonlinearities in the
boundary conditions contribute a further complication to hydrodynamic problems.
Such problems are usually governed by Laplace's equation (4.18) with free surface
conditions
(4.27)
ot oz OX ox oy oy
0<P = B_ g YJ _ J.... [(0<P)2 + (0<P)2 + (0<P )2] (4.28)
ot 2 ox oy OZ
in which 1] is the free surface elevation and B is a constant (the "Bernoulli
constant"). The first equation simply states that a fluid particle on the free surface
follows the free surface movement. The second equation is a consequence of having
a constant pressure on the free surface as expressed by Bernoulli's equation. Both
conditions are nonlinear. Further, there is no damping in the system so that any
input of energy from a numerical strategy or any entrapment of energy (e.g., in a
small, 2 Ax, wavelength or in a larger wavelength) leads to a divergent calculation.
In steady state calculations a stream function formulation is often superior to
that of the velocity potential. Using the stream function in (4.27) and (4.28) with
the time dervative set to zero, these equations become, respectively,
",=0 (4.29)
-1 (0",)2
- +g1]=B (4.30)
2 on
APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS 87
r- Zone of Yl
uncertainty--1
Thus the kinematic condition (4.27) has become linear. The dynamic condition
(4.28) is simplified and only contains the normal derivative, which is a natural
consequence of the calculation. Boundary derivatives in directions other than
normal are awkward to compute.
A typical problem is the two dimensional flow over a spillway (Fig. 4.4). The
solution region is terminated at the right and left boundaries where uniform
horizontal flow is assumed, olfl/on = 0, and IfI is linear along the boundary. Thus
the solution consists of finding the following ns + nf+ 1 unknowns (assuming a
known flow rate, q):
(b) Bernoulli's equation (4.30) at each free surface point (nf equations), and (c) the
specification of uniform flow through the left boundary, (olfl/on)s=(oy/on); (in
which the normal is with respect to the lower boundary and the free surface,
respectively).
Because Bernoulli's equation is nonlinear and because 1'/ is unknown, implying
nonlinearity in the boundary equations, the solution must be interative. After an
initial guess the solution can be improved by a standard Newton-Raphson
technique. Cheng et al. [39], however, found poor convergence in a zone behind
the spillway crest. Better convergence was achieved by creating a group of
"pseudo-nodes" on the free surface and satisfying Bernoulli's equation at the real
nodes plus the pseudo-nodes in a least squares sense. A similar result could also be
achieved by guaranteeing smoothness of the free surface through the use of cubic
spline interpolation functions on the free surface (Liggett and Salmon [40]).
88 APPLICATIONS OF BOUNDARY ELEMENT ME1HODS TO FLUID MECHANICS
(4.33)
in which p is the angle the free surface makes with horizontal. The subscript on the
derivative (oCP/ot)x indicates that this derivative is taken with the x-coordinate
held constant but the derivative follows the free surface vertically; that is, it is the
rate of change of cP on the free surface at a fixed horizontal point.
Equations (4.32) and (4.33) are suited for a time stepping scheme. Beginning
with an initial free surface at time kAt, (4.32) can be used to find the free surface
position at time (k+l)At and (4.33) gives the potential at time (k+ I)At. A
possible centered finite difference scheme is
+1 2(At) [(oCP)k + (1 -
r/' = r/' + cos pk + cos pk+1 (h a;; ( 1)
(oCP)k+l]
a;; (4.34)
(~~ L B- g ~ - ~ [ ( 00:
= 2( 1 f~
tan p + (4.36)
in whichfis the damping coefficient. Betts and Mohamad solved the problem of a
progressive wave in a channel with a numerical wave maker in one end. For a few
90 APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS
nodes at the end opposite to the wave maker they used equation (4.36) with/=F 0,
thus providing a wave absorber which reduced reflection. That, however, does not
solve the usual problem where the effects of a wave on a structure are studied. In
that case reflection from the structure would again reflect from the wave maker
and contaminate the solution. Better transmitting or absorbing boundary condi-
tions, which allow the outgoing wave to pass while not effecting the incoming
wave, are urgently needed.
In the small amplitude wave problem the free surface boundary conditions (4.32)
and (4.33) are linearized to give
017 orfJ
(4.37)
ot oz
orfJ
-=-g17 (4.38)
ot
in which z is the vertical coordinate. Moreover, these conditions are applied at the
equilibrium surface so that the point of application of the boundary conditons is
no longer a part of the solution.
Several different cases of linear wave problems are summarized herein.
Additional detail is found in the references below and in Liu and Liggett [41].
Assuming monochromatic waves of frequency OJ, the substitutions
17 = '(x, y) e- iwt rfJ = rp(x, y, z) e- iwt (4.39)
lead to
V 2rp = 0 (4.40)
with the free surface boundary condition
orp OJ2
-=-rp (4.41)
oz g
In the above equations i = 0 and the final solution results from taking the real
(or imaginary) part.
Two cases of two-dimensional flow are important. Vertical plane problems
assume a uniform behavior in the horizontal y-direction. In the case where rigid
surfaces are uniform in the y-direction
rp(x, y, z) = rp* (x, z) eik • y (4.42)
_ [i g A coshcoshk (zk +h h) e
II' (x, y, z) - - OJ
ikx
+ 11'* (x, y)
]
(4.45)
where HI and Ho are Hankle functions of the first kind of order one and zero.
For transient linear wave problems a time stepping technique along the lines of
that cited for nonlinear waves is easy to devise. Fortunately, the stability and
accuracy of such a linear scheme can be analyzed (Salmon et al. [54]). Three-
dimensional linear wave solutions are a solution of Laplace's equations under the
transient, free surface boundary conditions (Lennon et al. [55]).
Even in the linear case suitable radiation boundary conditions are difficult to
find and are need~d to limit the solution region. For shallow water such conditions
are given by Lennon et al. [55] but the shallow condition is inaccurate for deep
water waves. At a penalty in computing time the method of Smith [51] can be
used.
If all inertia terms are neglected in the Navier-Stokes equation and a stream
function formulation is used, the viscous flow statisfies the biharmonic equation
(4.48)
By use of the Rayleigh-Green identity the equation analogous to (4.2) (without the
volume integral) is (Jaswon and Symm [2], p. 271)
IfI ar 1 alfl a 2 2 ar ]
PIfl(P)=-J [ - - + - - + r - ( V 1fI)-(V 1fI)- dA (4.49)
r r2 an r an an an
in which P= 4n for pin Rand P= 2 n if p is on a smooth part of r.
Creeping flow inside a cylinder was computed by Mills [56] who used a more
elaborate Green's functions which was limitec;l to the cylinder. A general formula-
tion is given by Okabe [57).
92 APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS
4.9 Porous-Elasticity
by piecewise linear functions. (b) The amount of data which must be supplied is
minimized by boundary methods relieving the user of elaborate data input.
Similarly, but of less importance, it is easy for the user to choose the output that he
needs.
2. Special Situations. It appears to be easier to design boundary methods to handle
automatically singularities, moving boundaries and infinite boundaries.
3. Immediate Results. The variables of interest can be obtained from the method as
opposed to some domain methods where an additional numerical step must be
performed (e.g. in domain methods the potential flow velocity usually comes from
a numerical differentiation of the potential leading to a degradation of accuracy).
4. Efficiency.
5. Use on Mini- and Micro-Computers.
Not explicitly in the above list may be the biggest advantage of all - the ability
to solve huge problems. Even though computing is becoming less expensive and
computer storage more abundant, the shear size of many practical problems form
the greatest obstacle to their solutions. Several of the problems mentioned on the
previous pages fall into that category. Size may be a result of physical complexity
or of nonlinearities which require fine discretization, small time steps, and
repeated iterations. For that class of problem boundary methods may be the only
economically feasible technique.
References
I Trefftz, E., Uber die Kontraktion kreisfOrmiger Fliissigkeitsstrahlen. Z. Math. Phys., 64,
34, 1917
2 Jaswon, M A, Symm, G. T., Integral Equation Methods in Potential Theory and Elasto-
statics. Academic Press, 1977
3 Brebbia, C. A, The Boundary Element Method/or Engineers. John Wiley and Sons, 1978
4 Brebbia, C. A, Walker, S., Boundary Element Techniques in Engineering. Newnes-Butter-
worths 1980
5 Banerjee, P. K, Butterfield, R., Boundary Element Methods in Engineering Science.
McGraw-Hill 1981
6 Liggett, 1. A, Liu, P. L.-F., The Boundary Integral Equation Method/or Porous Media Flow.
George Allen Unwin, 1983
7 Crouch, S. L., Starfield, A M., Boundary Methods in Solid Mechanics. George Allen and
Unwin, 1983
8 Mukherjee, S., Boundary element methods in creep and fracture. Applied Science
Publishers, Ltd, Essex, 1982, Also Elsevier Science Publishing, Co, N.Y.
9 Green, G., An Essay on the Application 0/ Mathematical Analysis to the Theories of
Electricity and Magnetism, Nottingham 1828
10 Hunt, B. W., Numerical solutions of an integral equation for flow from a circular orifice.
Jour. Fluid Mech., 31, 361-377,1968
II Prager, w., Die Druckverteilung an Korpern in ebener Potentialstromung. Physik.
Zeitschr., 29, 865-869, 1928
12 Nystrom, E. 1., Uber die praktische Auflosung von linearen Integralgleichungen mit
Anwendungen auf Randwertaufgaben der Potentialtheorie. Soc. Sci. Fennica, Comment.
Physico-Math., 4, 15, I-52, 1928
13 Fredholm, 1, Sur une c1asse d'equations fonctionelles. Acta Math., 27,365-390, 1903
94 APPLICATIONS OF BOUNDARY ELEMENT ME1HODS TO FLUID MECHANICS
14 Lotz, I., Calculation of potential past airship bodies in yaw. NACA TM 675, 1932 [also
Ingenieur-Archiv, Vol. 11,1931]
15 Vandrey, F., A direct iteration method for the calculation of velocity distribution of bodies
of revolution and symmetrical profiles. Admiralty Research Lab. Rept. R IIG/HY 112/2,
1951
16 Weinstein, A, On axially symmetric flows. Quarterly of Applied Math., 5, No.4, 1948
17 Van Tuyl, A, On the axially symmetric flow around a new family of half bodies.
Quarterly of Applied Math., 7, No.4, 1950
18 Sadowsky, M. A, Sternberg, E., Elliptic integral representation of axially symmetric flows.
Quarterly of Applied Math. 8, No.2, 1950
19 Landweber, L., The axially symmetric potential flow about elongated bodies of revolu-
tion. David Taylor Model Basin Rep., 761,1951
20 Vandrey, F., On the calculation of the transverse potential flow past a body of revolution
with the aid of the method of Mrs. Flugge-Lotz. Astia AD-40089, 1951
21 Smith, A M. 0., Pierce, J., Exact solution of the Neumann problem. Calculation of plane
and axially symmetric flows about or within arbitrary boundaries. Douglas Aircraft
Company Report No. 26988, 1958 [Summary in Proc. of the 3rd Int. Congress of Applied
Math., Brown University, 1958]
22 Hess, J. L., Smith, A M. 0., Calculation of nonlifting potential flow about arbitrary three-
dimensional bodies. ES 40622, Douglas Aircraft Corp., Long Beach, Calif. 1962 (Also in
Jour. of Ship Research 8, No.2, Sept. 1964)
23 Davenport, F. J., Singularity solutions to general potential flow airfoil problem. D6-7207,
Boeing Airplane Co., Seattle, Wash. 1963
24 Jaswon, M. A, Integral equation methods in potential theory: I. Proc. Royal Soc. A, 275,
23-32,1963
25 Symm, G. T., Integral equation methods in potential theory; II. Proc. Royal Soc. A, 275,
33-46,1963
26 Chaplin, H. R., A method for numerical calculation of slip stream contraction of a
shrouded impulse disk in the static case with application to other axisymmetric potential
flow problems. David Taylor Model Basin Report No. 1857,1964
27 Gallagher, R. H., Finite Element Analysis Fundamentals, Prentice Hall 1975
28 Hess, J. L., Review of integral-equation techniques for solving potential-flow problems
with complicated boundaries. Innovative Numerical Analysis for the Applied Engineering
Sciences. University Press of Virginia, Charlottesville, 131 - 143, 1980
29 Hunt, B., The mathematical basis and numerical principles of the boundary integral
equation method for incompressible potential flow over 3-D aerodynamic configurations.
Numerical Methods in Applied Fluid Dynamics (B. Hunt, ed.), Academic Press, 49-135,
1980
30 Carey, G. F., Extension of boundary elements to lifting compressible aerodynamics. Finite
Element Flow Analysis (T. Kawai, ed.), Univ. of Tokyo Press, 939-943, 1982
31 Carey, G. F., Kim, S. W., Extension of boundary element method to iifting sub critical
flows. 19th Annual Meeting, Society of Engineering Science, University of Missouri-Rolla,
1982
32 Inamuro, T., Adachi, T., Sakata, H., A numerical analysis of unsteady separated flow by
discrete vortex model using boundary element method. Finite Element Flow AnalYSis
(T. Kawai, ed.), University of Tokyo Press, 931-938, 1982
33 White, J. W., Kline, S. J., A calculation method for incompressible axisymmetric flows,
including unseparated, fully separated, and free surface flows. Report MD-35, U.S. Air
Force Office of Scientific Research Mechanics Divison, Contract AF-F44620-74-C-0016;
Thermosciences Division, Department of Mechanical Engineering, Stanford University,
1975
34 Wu, J. c., Problems of general visous flow. In: Developments in Boundary Element
Methods - 2 (P. K Banerjee and R. P. Shaw, eds.). Applied Science Publishers, 69-109,
1982
35 Wu, J. c., Thompson, J. F., Numerical solutions of time-dependent incompressible
Navier-Stokes equations using an integro-differential formulation. Computers and Fluids
1,197-215,1973
APPLICATIONS OF BOUNDARY ELEMENT METHODS TO FLUID MECHANICS 95
62 Cleary, M. P., Fundamental solutions for fluid-saturated porous media and applications
to local rupture phenomena. Thesis submitted in partial fulfillment of requirements for
for the Ph.D. degree, Brown University, 1976
63 Cheng, A H-D., Liggett, 1. A, Boundary integral equation method for linear porous-
elasticity with applications to fracture propagation. Int. 1. for Numerical Methods in
Engineering, (In press), 1983
64 Cheng, A H-D., Liggett, 1. A, Boundary integral equation method for linear porous-
consolidation. Int. 1. for Numerical Methods in Engineering, (In press), 1983
65 Liggett, J. A, Hydrodynamic calculations using boundary elements. Finite Element Flow
Analysis (T. Kawai, ed.), Univ. of Tokyo Press, 889-896, 1982
ChapterS
by M C. Au and C. A. Brebbia
5.1 Introduction
Over the last decades the study of water waves diffraction and radiation problems,
especially in offshore structures has interested many researchers. The problem is
sometimes solved by using the finite element method as shown by Newton [1];
Berkhoff [2]; Bai [3]; Yue, Chen and Mei [4], and more recently Zienkiewicz and
Bettess [S]. Alternatively boundary integral equations have been used as illustrated
by the work of Garrison [6]; Hogben and Standing [7]; Eatock-Taylor [8] and
Isaacson [9].
The boundary element method [10] is a combination of classical boundary
integral techniques with finite element concepts and thus, combines some of the
advantages of the two methods. In this chapter boundary elements are used to
solve wave diffraction-radiation problems in a simple and general way. The
fundamental solution applied here are IIr or In (1/r) for three and two dimensional
problems rather than the more complex solutions applied by other authors [6], [7].
The advantage of this simple approach is that it becomes easier to write a general
program as those solutions are accurate and give stable results. The disadvantage
however is that more elements are usually required to solve a problem and because
of this taking advantage of the system symmetry becomes more important. Several
types of symmetry are presented here to indicate how to reduce computational
effort. The examples included in this chapter serve to illustrate the applicability
and potentialisies of the technique.
Governing Equations
For the wave structure interaction problem described in Fig. S.1, the fluid is
usually assumed to be incompressible and inviscid. Under these conditions the
fluid will remain irrotational and a velocity potential VI can be defined such that,
V=V' VI (S.l)
where V is the velocity vector. In this case the continuity equation can be written
as,
(S.2)
98 WATER WAVES ANALYSIS
surface
I
ouler boundary
The water surface condition for small amplitude waves can be expressed in function
of the potential as,
[PUt oUt
j)f2+g8;;=O (5.3)
where g is the acceleration of the gravity and n is the normal to the water surface.
A moving-body surface will produce a normal velocity given by,
oUt .
-=o·X (5.4)
on
where 0 is the unit normal vector to the surface and X is the velocity vector of the
structure.
Incident Wave
The time dependent surface elevation, ¢o, due to a monochromatic wave can be
written as
(5.5)
where ao is the surface incident wave amplitude; w is the wave frequency and Re
indicates the real part of the complex value between brackets. k, and k2 are the
wave number components in x and y directions.
The incident wave potential in the fluid can be expressed as
Uo = Re {uo e- iwt } (5.6)
where
Uo = igao coshk(z+h) e'"(k ,x +k ,Y)
- --
w cosh k h
with k, = k cos rx and k2 = k sin rx; where rx is the angle of incident wave with
respect to x-axis; h is the depth of the water; k is the wave number which can be
defined by the dispersion expression as
w2
ktanhkh=- (5.7)
g
Boundary Value Problem
The presence of the structure (Fig. 5.1) will disturb the incident wave and the
total potential of the fluid can be written in terms of the incident and diffracted
wave potential, i.e.
(5.8)
WATER WAVES ANALYSIS 99
with
U = Re (u e- iW1 )
where u is the diffracted wave potential which IS governed by the continuity
equation, i.e.
(5.9)
with the following boundary conditions:
(i) On the T] surface of the structure,
ou ouo
-+-=n'V onT] (5.10)
on on
(ii) On the horizontal sea bottom T 2 ,
ou
-=0 onT2 (5.11)
on
(iii) Sufficiently far from the body the diffracted outgoing wave can be approxi-
mated by the radiation condition:
ou
--iku=O onT3 (5.12)
on
The boundary T 3 should be vertical and at a distance away from the source of
disturbance.
(iv) The water surface condition (5.3) can be written as
ou w2
---u=O onT4 ( 5.13)
on g
T4 is the undisturbed water surface.
The above four different types of boundary conditions can be summarized by
the general expression
ou
-+pu-q=O onT (5.14)
on
where T = T] + T2 + T3 + T 4. P and q are parameters defined by the above
conditions (i) to (iv).
The boundary value problem defined by equations (5.9) and (5.14) can be
formulated in terms of a weighted residual expression [10] and solved approxi-
mately. Given an approximation for the u potential and a weighting or distribution
function u*, one can write,
Ci J ou* )
Ui + ( a;;- + p u* u dr = Ju* q dr (5.19)
where
Ci = t for the point 'j' on the smooth boundary r;
Ci = 1 if the point 'i' is in Q and
Ci = 0 if'i' is external to Q domain and boundary r.
r = L, re (5.20)
N
CiUi+L, S (0*)
oUn +pu* [<1>]dr{u~}=L, S u*[<1>]dr{q~}
N
(5.23)
e Te e re
\
/
I---- Radiation
boundary
~e sea bottom
Matrix Formulation
By considering that the point 'i' of equation (5.24) refers to every nodal point on r,
a complete matrix system of equations can be formed such as
[H] {u} = [G] {q} (5.25)
where [H] = [H] + [C]. Note that [C] is a diagonal matrix formed by the Ci
coefficients Of equation (5.24).
One can now apply the four different types of boundary conditions to equation
m-
(5.25) and obtain, after suitable rearrangement, the following system.
or
[A]{u} = {B}
where the matrices [Hj]ie and [Gj]ie are due to the integrals of the type
ou*
S -[(/>] dr and S u* [(/>] dr
~nG on ~nG
respectively.
The j in rj (j = 1,2, 3, 4) represent the four types of boundary conditions given
by (5.10) to (5.13) and the integrations are to be carried out over the element re,
usually using a system of local coordinates as in finite elements.
The {qd vector in equation (5.26) can be given by the incident wave and the
body surface velocities, i.e. at any point,
ql = qo + n' {"X} (5.27)
oUo. .
where qo = - IS gIven y t e mCI d ent f'Ie ld . I n matnx
b h . . ...lorm,
on
{qd = {qo} + [Q] {X] (5.28)
where
or
[Q] = [Gd [nx, ny, n" nyz , nzx, nxy] in 3-D
(nx, ny, nz) is the unit normal vector on r l
(nyz, nzx , nxy ) is {(x - xo), (y - Yo), (z - zo)} x {nx, ny, nz}
102 WAlER WAVES ANALYSIS
Wave Forces
Once the potentials around the structure are known, one can compute the pres-
sures by using the linearized version of Bernoulli's equation, i.e.
o
{P} = - e ot {U t} (5.32)
where e is the density of the fluid and Ut are the total potentials i.e. Ut = Uo + U.
For harmonic motion,
p = i e w (uo + u)
and (5.33)
P = Re {p e- iwt }
After the pressures are known the forces acting on the structure can be computed
by integration over the surface, i.e.
{f}=ie w S {n}Utdr
and 6xl r , 6xl (5.34)
{F} = Re [{f} e- iwt ]
Finally the water surface elevation can be obtained at any point. It is given by the
kinematic condition, which for harmonic motion gives,
~ = Re [ i w g (uo
. ]
+ u) e- 1wt (5.35)
In this section, several special cases will be studied, starting by the vertically
integrated structure, i.e. structures of constant horizontal sections and for which the
problem can be integrated with respect to depth; this effectively reduces the
problem to two dimensions. Other types of two dimensional structures are those
which can be studied on a vertical plane because they are sufficiently long in the
direction normal to that plane. Fully three dimensional structures are also
considered in detail, stressing the importance of their rotational symmetry.
WATER WAVES ANALYSIS 103
Cylinder
-
Incident wave
Water surface
Bottom
(i) Vertically Integrated Structures. This type of structure has a constant horizontal
section throughout the water depth as shown in Figure 5.3. The incident wave
potential - equation (5.6) - can be written as,
_ cosh k(z + h)
Uo = uo(x, y) cosh k h (5.36)
i gao ·(k k ) . .
where ilo(x, y) = - - - e' I H ,y. Provlded that the X - Y plane formulation for
w \
the part of the velocity potential independent of depth Z. The diffracted potential
can be similarly written as
_ cosh k(z + h)
u = u (x, y) cosh k h (5 .37)
Notice that (5.36) and (5.37) satisfy the boundary conditions on the water surface
and on the sea bottom. The problem can then be expressed in terms of the
diffracted potential il (x, y) and the modified boundary value problem for the x- y
plane becomes
in Q (5.38)
with boundary conditions,
oil + oilo = 0 on r
on on
and (5.39)
oil
--ikil=O on r
on
or in general oil
- + P il = q on r +r
on
By using the weighted residual technique one can write,
ou*
Su[(V2+k2)u*]dQ= S u ( -+pu* ) dr+ S u*qdr (5.41)
Q r+r' on r+r'
The fundamental solution corresponding to the above problem is given by the
solution of
(5.42)
and can be written as i
u* =-HO(k r) (5.43)
4
where r = Vex - x;)2 + (y - y;)2
Notice that solution (5.43) automatically satisfies the boundary condition on r', i.e.
ou*
lim S u ( -;)-+pu* ) dr=O (5.44)
r-+ 00 T' un
Now equation (5.41) can be written for every point 'i' on r as,
ou*
C; Uj+ S-;)- U dr = S u* q dr (5.45)
r un r
where r is the "internal" boundary usually defined by a clock-wise direction loop
over the structural section. The matrix form of equation (5.45) can now be
produced assuming that r is discretized into N elements. This gives
[H] {u} = [G] {q} (5.46)
where for a constant element,
H 'J.. = rS -on
j
l
o -4i Ho(k r) ] dr + b··IJ e-I
i
Gij = S - Ho (k r) dr
r 4
j
q. = - {-ou o}
J on J.
Notice that the elements of the vector {q} are equal to the incident potential
derivatives computed at the structural nodes.
Equation (5.46) can then be solved and the potential can be computed as
_ _) coshk(z+h)
Ut = ( Uo + U ------'.-----'- (5.47)
cosh k h
The total three dimensional forces on the vertical structure can be obtained by
I. I
integrating (5.47) in the vertical direction which gives
JlhIx =IQW
tanh k h
k
Jnx
Sln y (5.48)
/xy r nxy
WATER WAVES ANALYSIS 105
and
fyz}
{fzx l I kh fy}
= zO+T tanh 2 - fx
1{
where Zo is the centre of the moments /yz and fzx.
(ii) Structure with a Vertical Axis or Plane of Symmetry. A structure as the one
shown in Fig. 5.4 has a vertical axis of symmetry with two symmetrical parts r l
and r 2 • If the boundary conditions on the z-axis are not known, a full
discretization of boundary elements will be required in r I and r 2 This gives a
matrix equation of the form,
(5.49)
where the subindex I and 2 in {u} and {q} refers to values on r l and r 2
respectively. It can be seen by inspection that HII = H22 and H21 = HI2 and
similarly for the terms in [G).
11
I
where 's' refers to the symmetrical and 'a' to the antisymmetrical component.
Substituting (5.50) into (5.49) and rearranging the submatrices one obtains
l(HIl +0 H I2 ) 0
(HIl - H 12 )
] {us}
Ua =
l(G Il + G12
0
0
(Gil - G I2 )
] {qs}
qa
(5.51)
The full matrix equations presented in (5.49) can then be reduced to a diagonal
form as shown in (5.51). Therefore the solution for the symmetric and anti-
symmetric components can be obtained independently. Once they are computed
one can find the values of {ud and {U2}. However, if the wave forces are the only
result required, they can be obtained from {us} and {u a} as follows,
(5.52)
106 WA1ER WAVES ANALYSIS
Using this technique the computer time required to solve a two dimensional
problem can be reduced by approximately one third. For some three dimensional
structures, one or more planes of symmetry can be considered and much larger
savings in computer time achieved. Some of the examples presented in sec~ion 4
require only 15% of the solution time for the fully discretized problem. I.e. where
symmetry is not taken into consideration. For further details see reference [11].
(iii) Structures with Rotational Symmetry. Consider a domain with rotational
symmetry about the z-axis as shown in Figure 5.5. The total boundary surface r
can be represented by n
r=Lr l 1=1
(5.53)
If the boundary element mesh has the same type of symmetry (i.e. N elements on
every r l , 1= 1,2, ... , n) the matrix equivalent of equation (5.23) can now be
formed by considering the point i only on the r l segment. This gives
where the submatrices [hd and [gd are of dimension Nx N whose components are
given by
and
[hdij = S
rj
(ou*
on
+p u*) [cP] dr (5.58)
[gdij= S
rj
u* [cP] dr
where rj = C~ rJ. rj is the /h element in the fh symmetrical component,
1= 1, 2, ... , nand j = 1, 2, ... , N.
Alternatively [hJJij and [g Ik can also be defined by
[hdij = S
rJ
(ou*
on
+p u*) [cP] dr (5.59)
and
[gdij=
rJ
S u* [cP]dr
with the following transformation to the 'i' point
where the right hand side vector represents the coordinates of the i point on r l ,
1= 1,2, ... , n.
Considering now the full matrix for all positions of the 'j' point on the complete
boundary r, a matrix similar to (5.27) will be obtained, such that,
hI h2 h 3 ··· h n
hn hI h 2 · ··hn - I
[H]= hn - I h n hI··· h n - 2 (5.60)
h2 h3 h 4 ··· hI
Similarly for [G]. A transformation can now be applied to this type of matrix such
that [12]
{u} = [R] {u}
and (5.61 )
{q} = [R] {q}
where [R] is anN x n N matrix with submatrices
[Rim] = [I] ei (l-I)(m-I)8 for 1=1, ... ,n, m=I, ... ,n (5.62)
NxN NxN
108 WATER WAVES ANALYSIS
The summation defined in (5.66) is over n 2 tenns, i.e. I and m (both varying from 1
to n); each tenn [Him] in the equation needs to be considered once. If the order of
summation is altered in such a way that the direction of the summation will be in
the diagonal way, once can do this by replacing the subindex m by m' such that
{
m + 1- 1 if m' ~ n
(5.67)
m' = m + I - I - n if m' > n
The submatrix
(5.68)
for p= n. Therefore, using the submatrices defined in equation (5.60) one can
write,
[il kk ] = L [hm]ei(m-Ij(k-I)O (5.72)
NxN m NxN
WATER WAVES ANALYSIS 109
Hence the [ill matrix has been reduced to a new matrix with submatrices only on
the diagonal position; i.e.
illl 0 0 ... 0
0 il22 0 ... 0
[ill = 0 0 il33 •• • 0 (5.73)
0 0 0 ... ilnn
The same consideration can be applied to the matrix [G]. The final form of
equation (5.73) can be written as,
(5.74)
k = 1,2, ... , n and the dimension of the matrices are N x N. The unknown vectors
are given by
{ad = L [Rkjr l {Uj}
j
(5.75)
{qd = L [Rkjr {qj}
l
j
[Rr l =~[R*l
fJ
(5.76)
After obtaining the solution of equation (5.75) the wave forces can be calculated
from;
{f}=iQw S {n}Ut dr (5.77)
6xl r,6xl
[~ ~ ~J
nz nz
= (5.78)
nyz nyz
nzx nzy
nxy , nxy I
One can find the summation for the forces in (5.77) by considering an identity
similar to (5.70), i.e.
(5.79)
and the relation defined in (5.61). The summation becomes, after rearranging the
terms,
(5.80)
110 WA1ER WAVES ANALYSIS
where
0 (nx - i ny) I (-I(n~ nx+
+ i ny)
ny)
0 (inx+ny)
2nz 0 0
[N]=~
2
0 (n yz - i nzx) (nyz + i nzx )
0 (nxz + i nyz ) (nzx - i nyz )
2n xy 0 0
6x3
and {ud, {U2} and {un} are the solutions from the equation (74) for k = 1,2 and n
respectively. In such a way only three sets of equations are required to be solved.
The equations of motion for the structure can now be included in the analysis to
determine the movements of the system and how the forces change during motion.
Consider the following system
[M] {X} + [Cl {X} + [Kl {X} = {f} e- iOJ1 (5.81)
with {X} = {Xo} e- iOJ1• The mass matrix [M] is given by the structural mass matrix
[Mol and the added mass [Mal; i.e. [M] = [Mol + [Mal. The damping matrix is [C)
and together with [Mal can be obtained using (5.30) and (5.33) as
[Mal {X} + [C) {X} = - i Q OJ S {n} [Hr l [Ql dr {X} (5.82)
r,
[K] is the total stiffness matrix due to the contribution of the hydrostatic force and
those due to the mooring system. The time independent version of (5.81) is
{- OJ2 [M] - i OJ [C) + [K)} {Xo} = {f} (5.83)
whose solution is given by
{Xo} = [K - OJ2 M - i OJ C)-I {f} (5.84)
Several examples are considered here to illustrate the application of the above
theory and show the advantages of using boundary elements, expecially for those
cases where full advantage can be taken of symmetry. The examples cover a variety
of important problems and are
(i) Submerged half-cylinder,
(ii) Fully submerged cylinder,
WA1ER WAVES ANALYSIS 111
[ C] _-.:...,[C::-=-]_ (5.85)
3x~ = - e(n a2/2) w
Incident wave
& - (Water surface
~-----------5a------------~
ob
10- 1 4 6 8 10-' 1 4 6 8 1
ko-
Fig. 5.8. The wave forces on submerged half-cylinder by 50 elements. a Horizontal force. b
Vertical force
•
Radiation boundary
h Cylinder
Fully submerged
cylinder
Bottom
Bottom
Only 50 constants elements were used on half the total boundary as shown in the
figure. The results due to the swaying and heaving motion are given in Fig. 5.13
and compared against those by Newton [1].
(iv) Moored Floating Cylinder [14]. The moored floating cylinder of Fig. 5.14 has
been analysed assuming that mooring lines remain straight. This type of structure
has been used as a breakwater and is normally moored in the way showing in the
figure. Amplitude of motion computed using 60 constant elements are given in
Figure 5.15.
(v) Surface Waves on a Sloping Bottom. The case of surface wave amplitude was
studied in the example described in Figure 5.16. The bottom is sloping towards the
land and two distinctive parts were considered i.e. the straight and sloping
WATER WAVES ANALYSIS II3
2.5,-----------------,
- - - Naftzger's and Chakarbarti's solution
o 14 quadratic
elements
0.5
2.5,------------------,
2.0
1.5
0.5
ob
10- 1 I 4 6 8 10- 1 I 4 6 8 1 4
ka-
Fig. 5.11. Wave forces on the fully submerged cylinder obtained using 14 quadratic elements.
a Horizontal force. b Vertical force
Water su rface
Radiation
bound ary
I
J
t h
i
I
I
Bottom
~W=5a
~-----W-----~ h=6a
Fig. 5.12. A boundary element mesh for the half-floating cylinder
114 WATER WAVES ANALYSIS
1.25,-----------------, ,---------------------,
Heaving Swaying
1.00
- - Boundary element method
o Finite element method
by Newton (1]
11,---------------------------,
10
c
9
z
Water surface
-
Incident wave
Bottom
x
~----------------w----------------~
3.0.---------------------,
- - - Alliney's solution
° 70 constant elements
0:5
o 2 4 8 10 12 14 16 18 20
Distance from the land (x)
Fig. 5.17. The wave amplitude on the free water surface
portions. The horizontal sea bottom was defined as r 2 boundary and the sloping
part as rt. 70 elements were used to discretize the complete boundary and the
results for wave amplitude compared against those by Alliney [15] in Figure 5.17.
(vi) The Square Caisson. The square caisson of Fig. 5.18 has been studied by
Mogridge and Jamieson [16]. They computed the wave force by a circular cylinder
approach and compared them to the experimental results. In this chapter, the
theory described in section 4 (i) can be applied. 24 constant elements were used
and two extreme cases for IX angle of incidence were defined (IX = 0° and IX = 45°).
The results produced agreed with those of reference [16] as shown in Figures
5.19(a) and 5.19(b).
(vii) Elliptical Cylinder. Another case similar to (vi) is the elliptical cylinder of
Fig. 5.20, which was studied by Goda and Yoshimura [17] using an analytical
method. 16 quadratic elements have been applied here to discretize the boundary
and the angles of incidence of the wave was IX = 30 ° and IX = 60 0. Results are
given in Fig. 5.21 where they are shown to agree well with those of Goda and
Yoshimura.
116 WA1ER WAVES ANALYSIS
/x-
y
Incident wave
\
1111· x
a=6
hla=1013
2.5
.s 1.5
-"=
'"g; 1.0
.8r
.....x
0.5
0
0
3~'
15'
-15'
a>
~-3O'
~
-45' 0
- - 12 constant elements
a
o Mogridge's and
Jamieson's theory
-60'
• Mogridge's and 0
Jamieson's experiment
-75'
a b
-90'
0 0.2 0.4 0.6 0.8 1.00 0.2 0.4 0.6 0.8 1.0
kaht- kaln-
Fig. 5.19.
The magnitude and the phase of the horizontal force on a square caisson. a (X =0 0 ; b (X = 45°
WA1ER WAVES ANALYSIS 117
Incident wave
y
-fifff-
x
r
a=lO b=0.15a h=a
1.0.--------------~
~ 0.8 «= 30·
.....
~
..c:
.E 0.6
:§
-.J"
a
0
1.2
16 quadratic
- - elements
1.0
o Gada and
..c:::: Yoshimuro's theory
~
:;::
~ 0.8
g
c: «= 60·
"§
-: 0.6
<::>
~
-..;;
.......~
0.4
0.2
b
00 4 6 8 10
ka-
Fig. 5.21. Thefx force (a) and the!;. force (b) on elliptical cylinder by quadratic element
118 WATER WAVES ANALYSIS
(viii) Submerged Storage Tank and Tower. The geometry of this three dimensional
structure is given in Fig. 5.22 and in Fig. 5.23 the discretization shows how the
planes of symmetry at 60 0 have been used. The 56 constant triangular elements
applied on the one-sixth of the total boundary were used to define the structural
surface, the sea bottom, water surface and radiation boundary. Results obtained
using boundary elements are compared in Fig. 5.24 with finite element solution
due to Zienkiewicz and Bettess [5].
Tower
Storage tonk Radiaton Structural
surface surface
h
Fig. 5.23. The discretization for the submerged storage tank and tower with the relevant
boundary conditions
I
~ 0.08
f
~0.16
I
'1:: 0.04
~
8: ~ ~
-..!:' ....:-' Et
0.04 0.08
o 2 4 o 2 4 o 2 4
a kh- b kh- c kh-
Fig. 5.24. The wave forces on the submerged storage tank and tower. a Horizontal force. b
Vertical force. c Bending moment about the bottom
WATER WAVES ANALYSIS 119
(ix) Floating Hemisphere. This example has been solved by Garrison [6] using a
source distribution technique with a complex fundamental solution. In the present
case only the geometry of a vertical section was required and discretization was
applied as in Fig. 5.26, making use of rotational symmetry. The three dimen-
sional boundary element mesh for one symmetrical component covers 36 0 and can
be easily obtained by hand or by using an automatic mesh generating code. The
results obtained using 30 constant triangular elements are given in Fig. 5.27 and
made non-dimensional by the factor (lie gao a 2).
(x) Surface Elevation Around a Cylinder. The surface wave elevation around the
vertical cylinder in Fig. 5.28 was considered next. In this example, the values of
the potential on the water surface are required and it would be difficult to work
with a complex fundamental solution as in [6]. By using the theory presented in
h
~ ~
I
h
j
WAYJ.&-'<S,?&WY-<'Y,07$Y/'&/,M'N/,U/:l;J.&J.&/,&l/,(C
h=3a a=lO W I W=5a
Fig. 5.25. The floating hemisphere at the Fig. 5.26. A vertical section for the floating
free surface hemisphere
3.0.----------------,
- - Garrison's solution
2.5 - -- Boundary element
method
2.0
1.0
c: 2a
::1
h
x
h~a ~1 k~2
1.5,---------------0---------,
1.0
0.5
-1.0
--Theory
-1.5 0 Boundary element method
(rotational symmetry)
section 4 (iii), the water surface elevation was obtained using equation (5.34) and
40 elements were used due to symmetry. The results for the real and imaginary
part of the wave elevation are compared in Fig. 5.29 with the theoretical values.
The surface elevation was calculated using (5.37) i.e.
Free surface
Tower
Bose
Fig. 5.30. The boundary element mesh for the Condeep type structure
a
0.4 0.4 0.20 a
a
'&0.3 a
Q"
"-
'<-!"0.2
0.1 0.05
5.7 Conclusions
This chapter shows how the boundary element method can be applied to solve
wave diffraction-radiation type problems. By using the fundamental solution
corresponding to Laplace's equation results such as pressure, wave forces, hydro-
dynamic coefficients, water surface elevation and motion response of structures can
be easily obtained.
Two dimensional problems can generally be solved without great usage of
computer time. However, computer time becomes crucial when solving three
dimensional cases. For such cases the codes can be made more efficient and easy to
use (i.e. less data is required) by exploiting the symmetry of the structure. The
chapter shows how this can be done in practice.
122 WA1ER WAVES ANALYSIS
Several examples are considered at the end of the chapter to \ illustrate the
applications of the theory presented and to demonstrate some of the special
advantages of using boundary elements for these cases. The examples cover a
variety of important two and three dimensional problems, with and without
symmetry.
References
I Newton, R E., Finite element analysis of two dimensional added mass and damping.
Ch. 11, In: Finite Elements in Fluids - Vol. 1 ed. by R H. Gallagher et al. A Wiley-Inter-
science publication, 219-232,1975
2 Berkhoff, 1. C. W., Computation of combined refraction-diffraction. Ch. 24, 13th Conf. on
Coastal Eng., 1,471-490, 1972
3 Bai, K 1., Diffraction of oblique waves by an infinite cylinder. 1. Fluid Mech., 68,
513-535, 1975
4 Yu, D. K P., Chen, H. S., Mei, C. c., A hybrid element method for diffraction of water-
waves by three-dimensional bodies. Int. 1. Num. Meth. Eng., 12,245-266,1978
5 Zienkiewicz, O. c., Bettess, P., Fluid structura dynamic interaction and wave forces - An
introduction to numerical treatment. Int. 1. Num. Meth. Eng., 13,1-16, 1978
6 Garrison, C. 1., Hydrodynamic loading of large offshore structures: three dimensional
source distribution methods. Ch. 3, Numerical Methods in Offshore Engineering, ed. by
O. C. Zienkiewicz et aI., A Wiley-Interscience Publication, 1978
7 Hogben, N., Standing, R C., Wave loads on large bodies. Int. Symposium on the
dynamics of marine vehicles and structures in waves, Inst. Mech. Eng., 258-277, Lon-
don 1974
8 Eatock-Taylor, R, Generalised Hydrodynamic Forces on Vibrating Offshore Structures by
Wave Diffraction Technique. Offshore Structures Engineering edited by F. L. L. B. Carneiro
et al., Pentech Press, London, 249-274,1977
9 Isaacson, M de St. Q., Vertical cylinder of arbitrary section in waves. J. Waterway, Port,
Coastal and Ocean Div., ASCE, 309-324,1978
10 Brebbia, C. A, The Boundary Element Method/or Engineers. Pentech Press, London 1978
11 Au, M c., Brebbia, C. A, Computation of Wave Forces on three dimensional offshore
structures. Proc. of 4th Int. Seminar on Boundary Elements, Southampton and Boundary
Element Methods in Engineering. Springer-Verlag Berlin, Heidelberg, New York 1982
12 Zheng Xiaozhong, Bao Gang, Sun Shuxun, Application of Group theory to vibrational
analysis of shell structures with space rotation symmetry. Proc. Int. Conf. of Finite Ele-
ment Methods, 1982, Shanghai, China, edited by He Guangquian and Y. K Cheung
13 Naftzger, R A, Chakrabarti, S. K, Scatterin of waves by two-dimensional circular
obstacles in finite water depths. 1. Ship Research, 23, No. I March, 32-42, 1979
14 Ijima, T., Chou, C. R, Yoshida, A, Method of analysis for two-dimensional water wave
problems. Pap. 156, Proc. 15th Coastal Eng. Conf., ASCE, 2717-2736
15 Alliney, S., A numerical study of water waves on sloping beaches. Appl. Math. Modelling,
5, Oct. 321-328, 1981
16 Mogridge, G. R, Jamieson, W. W., Wave forces on square caissons. Pap. 133, Proc. 15th
Coastal Eng. wlonf., ASCE, 2271-2289
17 Gada, Y., Yoshimura, T., Wave forces on a vessel tied at offshore dolphins. Pap. 96, Proc.
13th Coastal Eng. Conf. Vol. IV, 1723-1742, 1972
18 Garrison, C. 1., TflJrum, A, Iversen, c., Leivseth, S., Ebbesmeger, C. C., Wave forces on
large volume structures - a comparison between theory and model tests. Proc. 7th Off-
shore Technology Conf., Paper No. 2137, Houston 1974
Chapter 6
6.1 Introduction
to heat flux, there is a difficulty in nodal assignment. Also, a difficulty arises when
using subdomains at points where several elements meet [4].
These difficulties arise when using continuous elements because freedom nodes
must be taken at the singular locations. If complete interelement continuity is
abandoned then freedom nodes may be moved inside the boundary element so that
they do not occur at the previously cited locations. The above mentioned problems
then disappear.
In this chapter suitable discontinuous elements are presented for two and three
dimensional applications and case studies are given to demonstrate their efficacy.
Additionally, in the case of three dimensional problems abandonment of inter-
element continuity allows freely graded, mismatched, meshes to be used advanta-
geously.
Boundary Equations
The Laplacian problem is typical of elliptic problems and is used here as an
example of. the development of boundary integral equations and of boundary
elements. Suppose a problem is defined on a domain Q with boundary r. Then the
problem is stated by the requirement that the governing homogeneous field
equations be satisfied in Q together with appropriate boundary conditions viz:
V 2 u (x) = 0, X E Q (6.1)
u(x) =u(x),
(6.2)
au au
- (x) = - - (x) X E r2
an an
Here r] + r 2 = r, the derivative is with respect to the normal, and u and au are
the given boundary values. Let w satisfy the governing equation in Q; an
(6.3)
Then
(6.4)
and after integrating equation (6.4) by parts twice and taking account of equation
(6.3) it follows that
ow au
- Su - dr + S- w dr = 0 (6.5)
r an r an
where the discrete term is the Cauchy prinipal value of the singular integral and c
is a known constant dependent on the dimensionality of the problem and the
regularity, or otherwise, of the surface at x'. On taking account of the known
boundary values (eqn. (6.2)), eqn. (6.8) becomes
ou* ou ou*
- S u-dr+ S-u*dr-Cu(x')= S i i - d r - S -u*dr (6.9)
au
r 2 on r , on r, on r on 1
where: the right hand side integrals are known, the discrete term is known if x E r1
and not otherwise, and the left hand side integrals involve the unknown u, °ou . On
Discretization
In modelling a planar problem the boundary is partitioned into a series of
segments which are either straight or curved. The straight segments are most easily
defined but usually at the cost of poor geometric representation of the problem
geometry. Constant or polynomial shape functions for the freedom pairs u, ~ (or
on
displacement and traction in a stress problem) are easily defined (Fig. 6.1 a).
!;
/C> /c
I
/"!;~+1
2" !;~ 0
x
3~
Geometric nodes
Freedom nodes
1 " C~-l 1" !; ~ -1
(i) Constant (ii) Linear
a b
Fig. 6.1 a Straight line element. b Quadratic geometric element
MI CO = t ((( - 1)
M 2 CO = 1 _ (2 (6.11)
M\() = t ((( + 1)
An ascending hierarchy of shape functions may be used [10] for the freedom
functions. Here a linear variation of the field quantities u, ~ is assumed. Then on
on
taking freedom nodes at the element ends (( = ± 1) the assumed forms of the
freedom functions are
u CO = FI UI + F2 U2
an CO
OU
= FI an + an
OUI
F2
OU2
(6.12)
FI = t (l - 0; F2 = t (l + 0 (6.13)
INTERELEMENT CONTINUllY IN THE BOUNDARY ELEMENT MElHOD 127
Suppose that the boundary is discretized into N elements, then the discretized
form of the singular boundary integral equation (Eqn. 6.8) is
N ou* ouW
ci ui + L, S uV)-dr = L, S - u* dr (6.14)
j=1 rj on j=1 rj on
ou* . . ou v )
where u*, -0- is the fundamental solution with singular point at x'; u U), - - are
n ~
the interpolation field functions of the ph element r j .
As previously mentioned it is not necessary to locate the singular point of the
kernel function u* on the boundary. If instead it is located at some exterior point
the discretized form of the boundary integral equation is
N 0* N oU)
L, S u v)_u_ dr = L, S _u_ u* dr (6.15)
j=1 r j on j=1 r on j
Now the integrals are no longer singular and higher order kernel functions may be'
used.
On expressing the field functions explicitly, as in eqn. (6.12), the regular
boundary integral equation becomes
l J
r
N N - -
OUI
j~ J; [FI F2l a;; dr uJ j~1 J; [FI F2l u* dr
ou* UI on
= OU2 (6.16)
. on
or eqUIvalently
(6.17)
where
ou* ou*
hi,j-I = S F I -
on dr; hi,j= S F 2 - 0- d r
r;_l r;(6.18)n
gi,j-I = S FI u* dr; = S F2 u* dr,
gi,j
r;_l r j
a, global numbering system has been applied to the nodal freedoms and the index
'i' labels the kernel function. In equation (6.17) N of the freedom variables are
prescribed by the boundary conditions, one at each node, so that the equation is a
non-homogeneous linear algebraic equation in the remaining N unknowns. On
employing N linearly independent kernel functions a determinate linear algebraic
system is obtained for the unknown paramters.
The procedure outlined appropriate to continuous elements is basically very
simple but it encounters difficulties as discussed in the next section.
Modelling Difficulties
When continuous elements are employed freedom nodes are located at the element
edge. Usually this poses no special problem but the element edge may coincide
with a discontinuity in the problem definition resulting in modelling ambiguities.
Thus in a two dimensional problem at a corner the normal to the boundary is
not defined; nevertheless in a Laplacian problem the freedoms are the field
128 INTERELEMENT CONTINUITY IN TIIE BOUNDARY ELEMENT ME1HOD
V
"
'-pointof
excessive
3 constraint
Freedom Interface
constraint
(I)
[} au _ au
an-Tn
a b c
Fig. 6.2. a Non-smooth surface. b Freedom constraint at an interface where boundary
condition type changes. c Freedom constraint at a point where several surfaces meet
function u and its normal derivative :: . Clearly the normal derivative is not
defined at the comer. Other elliptic problems exhibit the same feature since the
boundary equations result from integrations by parts, for instance in an elastostatic
problem the surface traction is not defined at a comer. In practice two approaches
·
have been adopted. First the problem has b een Ignored . au
In that a;; has
b een
treated as a scalar variable at the comer with resulting erroneous results (an
example of this appears in ref. [11]). Second the mesh has been rendered
discontinuous at the comer by not joining adjacent elements and leaving the
freedom nodes close together on either side Of the geometric singularity [1]. Good
approximate numerical results can then be obtained at the cost of incomplete
modelling - the comer region is ignored. But in the singular method equations are
generated using fundamental solutions with singularities at both adjacent nodes in
tum which is numerically objectionable since strong linear dependence is induced
in the algebraic system because of the effective linear dependence of kernel
functions. In three dimensional problems the same feature arises at a boundary
edge or cusp (Fig. 6.2a).
Another problem arises at the interface between F) and F 2 , that is where the
type of boundary condition changes. A freedom node on this interface is such that
for one element is known and for the other :: is known (see Fig. 6.2b).
Yet again, when the numerical model is derived using subdomains [4] a similar
problem arises where several subdomains meet or where two sub domains meet the
exterior boundary (see Fig.6.2c). The latter problems are discussed by Lachat,
Chaudonneret and Wardle et aI., [4], [12], [13].
Discontinuous Elements
Here, the geometry of the problem is defined, as usual, using straight line or
curved elements. Again, a hierarchy of linear or other higher order variation of the
field function and its flux-like variable may be taken. However, these freedoms are
not defined at the geometric nodes but at "freedom nodes" which can, in principle,
be freely chosen anywhere over the element. In the present assignment these
freedom nodes are chosen so that in a mesh of congruent elements these nodes are
equally spaced. For simplicity, the formulation of the planar problem using a
discontinuous linear element is presented here, followed by a brief description of
other higher order elements.
(6.19)
The algebraic expression of the boundary integral equation [Eq. (6.15)] for kernel
function 'i' now reads
.f.
J~1
(h;,2j-l U2j-l + h;,2j U2j) = f.
J~1
{g;,2j -l ( ~u).
un 2J-l
+ g;,2j ( ~u)
n
.}
2J
(6.20)
~
!
IX 1
~=+3/4
)/ ~=+1/4
where
au* au*
h i,2j-1 =J FI-a-dr;
r n
h i,2j = JF2 -a-
r n
dr
j j
(6.21)
gi,2j-1 = J FI u* dr; gi,2j = JF2 u* dr
r j rj
and the summationj is over each element with two nodes per element. Of the 4N
variables Uj, (~~) , 2N are fixed by boundary conditions, one at each freedom
node. On taking 2N linearly independent Kernel functions a determinate non-
homogeneous linear algebraic system is obtained.
(b) Discontinuous Quadratic Elment (Fig. 6.3 b). On taking equispaced freedom
nodes as before, interpolation functions are:
FI = ((9/8 ( - 3/4)
F2 = (1 - 3/2 0 (1 + 3/2 0 (6.22)
F3 = (9/8 (+ 3/4)
(c) Discontinuous Cubic Element (Fig.6.3c). The interpolation functions for the
four freedoms are:
FI = «( - «(
114) + 114) (1 - 4/30
F2 =(4 (- 1) «( - 3/4) «( + 3/4)
(6.23)
F3 =(-4 (- 1) «( - 3/4) «( + 3/4)
F4=«( - 114) «( + 114) (1 + 4/30
The implementation of these elements and generation of a determinate algebraic
system follows as for the linear element.
Type 1
Type 2
a b c
Fig. 6.4. a Partially discontinuous linear element. b Partially discontinuous quadratic
element. c Partially discontinuous cubic element
(b) Partially Discontinuous Quadratic Element. Again there are two variants,
Fig. 6.4 b, with interpolation functions:
Type 1 Type 2
(c) Partially Discontinuous Cubic Element. Appropriate to Fig. 6.4c, the interpola-
tion functions are:
Type 1 Type 2
aT =0 T=0
T=300 an 6
6
a x
-24.1
-51.8
-53.2
T~ aT
an 124.9
75
24.9
2251 175.1
b 300 c 275.1 .
-50.13 -50.0
-49.9
o Geometric nodes
200 200.01 x Freedom nodes
d 298.1 e 275
Fig. 6.5. a Numerical model for the prism problem. b Solution obtained using continuous
linear elements. c Solution obtained using discontinuous linear elements. d Solution obtained
using continuous linear elements with double nodes at corners. e Solution obtained using
partially discontinuous elements
Applications
Three two-dimensional problems, thermal, fluid flow and stress, are considered
which illustrate the relative merits of continuous and discontinuous elements. For
each problem the solution has been obtained using, (i) continuous elements, (ii)
continuous elements with double nodes at geometric singularities, (iii) dis-
continuous elements and (iv) partially discontinuous elements at corners and load
discontinuities with continuous elements elsewhere. Two noded linear elements are
used throughout. In each problem, as in the original references, cases (i) and (ii)
are examined in conjunction with the conventional (singular) boundary element
method while in cases (iii) and (iv) the regular boundary element method is used.
All these problems subsequently have been analysed for all four cases using both
INTERELEMENT CONTINUllY IN THE BOUNDARY ELEMENT METHOD 133
methods with the result that both methods yield sensibly the same results in each
case, the variations being due solely to case considered.
(a) Steady Heat Conduction through a Rectangular Prism. In this problem two
opposite faces of the prism are insulated while the other faces are held at different
constant temperatures giving a linear temperature variation between the non-
insulated faces. While this problem is very simple it has, nevertheless, been of
considerable use in tests [1], [11]. The numerical model employed, Fig. 6.5a, has
three equal linear elements on each face. Computed results for the four cases
considered are given in Figs. 6.5b-e.
(b) Steady Inviscid Laminar Flow in a Channel Past a Disc. A finite channel is
considered with the disc located symmetrically in the channel with its normal along
the channel. Only a quarter of the domain need be analysed on account of the two
symmetry planes present. A graded 32 linear element mesh, Fig. 6.6 a, is taken
together with a stream function, !fI, representation. Since the flow velocity is
infinite at the tip of the disc, this problem is singular [14], and so it presents a
severe numerical test. Computed values of the stream velocities at the channel
wall, B C, and over the disc face, DO, are given in Figs. 6.6 band 6.6 c respectively.
(c) Stress Concentration in a Rectangular Plate Having a Circular Hole. The
rectangular plate, 150 mm x 110 mm, is perforated at its centre by a circular hole of
radius 15 mm and is subjected to uniform tensile loading. Again, the symmetry
planes permit an analysis on the quarter geometry which is modelled using 36
graded linear elements, Fig. 6.7 a. The stress concentration factor for this problem
is known [15]. Computed values of the normal stress along the ligament ED are
given in Fig. 6.7b. The implied stress concentration factor is accurate for cases (ii),
(iii) and (iv) and is 7% low for case (i).
Discussion
The poor performance of the continuous elements at corners is abundantly evident
in the problems considered here. Nevertheless, the effect is localized, with no
general degradation being apparent so that there is no severe detriment to the
interior solutions.
In the heat conduction problem, the case of continuous elements with double
nodes at a singularity, case (ii), presents about 1% error in temperature at the
corners, A, B, and 0.2% error away therefrom. This is accompanied by 0.1 % error in
flux on the faces AD, CB. The temperature values obtained with discontinuous
elements, case (iii), are similar to case (ii), while the flux values show some 4%
error at the corners with 1% error away therefrom. A similar result has been
observed with constant elements [1]. The partially discontinuous elements, case
(iv), perform best with results within 0.1 % error over the whole boundary.
In the fluid flow problem, cases (ii), (iii), (iv) give the same results within 1%
on the channel wall and 2% on the disc face away from its edge. At the edge of the
disc, where the solution is singular, cases (iii) and (iv) give the same result with
cases (i) and (ii) trailing by more than 25%. In order to achieve results with case
(ii) or comparable quality to those of cases (iii) and (iv), the mesh density on the"
singular edge needs to be doubled.
134 INTERELEMENT CONTINUIlY IN THE BOUNDARY ELEMENT METHOD
I
A - - Continuous elements
B ----- Discontinuous elements
C .-6-.. . . . Partially discontinuous elements
Continuous elements with double nodes at corners
::tOA-A-I>--I>
'i;.O.5
co
~---O--~r--O--~~~~~~~~~~~-O~
b B c
Fig. 6.6. a Boundary element discretization of flow problem. b Stream velocities at the
channel wall, Be. c Stream velocities along the disc
1 elements
0
-0
L __
'"c
I
u
0 c EII>
5
Ux=0
ty =0 tx =1
ty =0 E
a tx =O;uy=O B b tx-
Fig. 6.7. a Boundary element discretization of plate problem. b Normal stress distribution
along the ligament "DE" of the plate
INTERELEMENT CONTINUITY IN THE BOUNDARY ELEMENT METHOD 135
The stress concentration factor given with cases (ii), (iii) and (iv) agrees with
Peterson's value within 0.3% while in case (i) the result is 7% low.
It is clearly evident that, while the discontinuous element models have around
twice the freedoms of the others, there is no consequent improvement in quality of
result. This is because, despite the increase in degrees of freedom, the mesh is no
more complete than a comparable mesh of continuous elements. Thus, discon-
tinuities should only be introduced where they are required for modelling
purposes.
Other case studies with singular problems [11], [16], have revealed that if the
mesh is sufficiently coarse in the neighbourhood of a singularity, then a spurious
oscillation can occur in the surface solution using continuous or discontinuous
elements. However, this disappears on mesh refinement and is not seen in the flow
problem discussed here.
where
MI (0 = 114(1 + (1)(1 + (2)«(1 + (2 - 1)
(6.28)
M 5(O = 112(1- (r)(1 + (2)
and the remaining shape functions follow by induction.
As usual, having specified the geometric element, there remains a wide latitude
of choice of freedom functions. Here attention is confined to 4, 8 and 12 freedom
136 INTERELEMENT CONTINUI1Y IN TIIE BOUNDARY ELEMENT MElHOD
~2
o Geometric nodes
3
Fig. 6.S. 8-node isoparametric quadratic lateral element
node variants. In each case the field quantities u, ~~ are given by:
n
u«) = L, Fj«) uj
j=1
(6.29)
iJu n. iJu j
-«)= L, £1«)-
an j=1 an
where the element has n freedom nodes, Fj (0 are the appropriate interpolation
functions and (= «(I, (2); (i E [- 1, + 1]; [17], [18].
Discontinuous Elements
(a) Linear Elements. Here, (Fig. 6.9a), the freedom nodes are taken at the four
points ( = (± i-, ± i-) and the interpolation functions are:
Fj «) = (i- ± (I)(i- ± (2) (6.30)
(b) Quadratic Element. The eight freedom nodes are located as shown in Fig. 6.9b
and typical interpolation functions are:
FI (0 = 27/32 (2/3 + (I) (2/3 + (2) «(I + (2 - 2/3)
(6.31)
F\O = 27/16 (4/9 - (1) (2/3 + (2)
the rest following by induction.
b c ~--<o---r--
(c) Cubic Element. Again, the twelve freedom nodes are located as in Fig.6.9c,
typical interpolation functions are:
FI (Q = 1/9(3/4 + (d (3/4 + (2) (8 «(r + (~) - 5)
(6.32)
P(Q = 3/2(1/4 + (1)(1- 16/9 (f)(3/4+ (2),
the rest following by induction:
b) Two Adjacent Edge Discontinuities. The freedom nodes are given in Fig. 6.lOb
and the interpolation functions are:
F1 eo = %(1 + (1) (1 + (2)
F2 eo = %(-i. - (d (1 + (2)
F3 «() = %(1 - (1) (1 - (2) (6.34)
reo = %(1 + (1)(1 - (2)
c) Two Non-adjacent Edge Discontinuities. The freedom nodes are given in Fig. 6.10c
and the interpolation functions are:
(6.35)
d) Three Edge Discontinuities. The freedom nodes are given in Fig. 6.IOd and the
interpolation functions are:
F1eo = t (1 + (d (1 + (2) (6.36)
reo = t (1+ (1)(1- (2)
with F2, F3 following by induction.
Applications
As yet, there is little published information regarding the use of discontinuous, or
partially discontinuous, elements in three dimensions. Discontinuous elements and
their applications, for the regular boundary element method have been dis-
cussed by Patterson and El-Sebai [17], and it has been reported by Danson [19],
that discontinuous elements are available in a commercial package, using the
singular boundary element method, but without case study material. In reference
[17] 4, 8, and 12 freedom node discontinuous elements were presented, in
conjugation with a 8 node isoparametric geometric element. The accompanying
case studies showed that the elements could be successfully applied, including
cases where subdomains and freely graded meshes were employed. Two of the
cited examples are considered here and the results are extended to include the use
of partially discontinuous elements.
300 ,-----------,
t 200
y = 1.25
"- 100
Z = 0.5
c X-
300 300
- Known solution
4
b Z- d y-
Fig. 6.11. a Bar problem - geometry and boundary conditions. b Temperature distribution
along AB. c Temperature distribution along CD. d Temperature distribution along EF
aTlan
~~~----I·I T=300
____________ .L_
I /
I I
1/
_ . _ . _ . - -W--.
\
\
-z
\
" T=O
a
[
T
" 1"",-
1\ I I
II -------------+-{
W-- _____________D+-\ B
Fig. 6.12. a Cylinder problem. b
V I \ Cylinder; problem-discretization.
/ 1/ c Computed surface values along
b
A ABCD
140 INTERELEMENT CONTINUITY IN 'IRE BOUNDARY ELEMENT MElHOD
half face is held at temperature 300 and the remaining surface of the cylinder is
held at temperature O. The solution of this problem clearly has an axial plane of
symmetry, so that only half the cylinder need be modelled, Figure 6.12 a. Further-
more, the solution is fully three dimensional and is singular at all edges of the
heated surface. A fourier series solution is available [20]. The element mesh em-
ploys twenty mismatched elements, Figure 6.12 b. Computed surface values on the
end-face are given in Figure 6.12c. The internal solution, on several internal lines, was
compared with the series solution and again there was less than 0.1 % discrepancy.
Both of these problems have been re-analysed using partially discontinuous
elements, with eight geometric and eight freedom nodes, giving results of similar
quality.
In view of the manifest incompleteness of the coarse meshes employed, in their
ability to describe the singular edge heat fluxes, the quality of the interior results is
remarkably good for both problems. The bar problem results show that both the
discontinuous and the partially discontinuous elements are well suited to sub-
domain modelling. Also, the cylinder problem results similarly validate the use of
freely graded meshes.
6.5 Oosure
References
1 Brebbia, C. A, The Boundary Element Methodfor Engineers. Pentech Press, London 1978
2 Watson, J. 0., Advanced implementation of the boundary element method for two and
three-dimensional elastostatics, In: Developments in Boundary Element Methods - 1. Eds.
Banerjee, P. K & Butterfield, R., Applied Science Publishers, London, Chapt. 3, pp.
31-63, 1979
INTERELEMENT CONTINUIlY IN THE BOUNDARY ELEMENT METHOD 141
3 Patterson, C., Sufficient conditions for convergence in the finite element method for any
solution of finite energy. In: Mathematics of Finite Elements & Applications. Ed. White-
man,1. R., Academic Press, London, pp. 213-224,1973
4 Lachat, 1. c., A Further Development of the Boundary Integral Technique for Elastostatics.
Ph.D. Thesis, University of Southampton 1975
5 Patterson, c., Sheikh, M. A, A regular boundary elementmethod for fluid flow. In: Int.
1. Num. Methods in Fluids, 2, pp. 239 - 251, 1982
6 Patterson, c., Sheikh, M. A, Regular boundary integral equations for stress analysis. In:
Boundary Element Methods. Ed. Brebbia, C. A, Springer-Verlag, Berlin, pp. 85 - 104, 1981
7 Patterson, c., Sheikh, M A, A regular boundary element method for coupled sub-
domains. In: Numerical Methods for Coupled Problems. Eds. Lewis, R. W., Bettess, P., and
Hinton, E., Pineridge Press, Swansea, pp. 115-126, 1981
8 Patterson, c., Sheikh, M. A, A higher order boundary element method for fluid flow. In:
Finite Element Flow Analysis. Ed. Kawai, T., Tokyo University Press, Tokyo, pp. 907-913,
1982
9 Ergatoudis, 1., Isoparametric Finite Elements in Two and Three-Dimensional Analysis.
Ph.D. Thesis, University of Wales, Swansea 1968
10 Brebbia, C. A, Walker S., Boundary Element Techniques in Engineering. Newnes-Butter-
worths, London 1980
11 Patterson, c., Sheikh, M. A, Discontinuous boundary elements for heat conduction. In:
Numerical Methods in Thermal Problems - II. Eds. Lewis, R. W., Morgan, K, and
Schrefler, B. A, Pineridge Press, Swansea, pp. 25 - 34, 1981
12 Chaudonneret, M., On the discontinuity of the stress vector in the boundary integral
equation method for elastic analysis. In: Recent Advances in Boundary Element Methods.
Ed. Brebbia, C. A, Pentech Press, London, pp. 185-194, 1978
13 Wardle, L. 1., Crotty, 1. M, Two-dimensional boundary integral equation analysis for
non-homogeneous mining applications. In: Recent Advances in Boundary Element Methods.
Ed. Brebbia, C. A, Pentech Press, London, pp. 233 - 251, 1978
14 Milne-Thompson, L. M., Theoretical Hydrodynamics. Macmillan, London 1968
15 Peterson, R. E., Stress Concentration Factors. John Wiley & Sons, New York 1974
16 Patterson, c., Sheikh, M. A, Non-conforming boundary elements for stress analysis. In:
Boundary Element Methods. Ed. Brebbia, C. A, Springer-Verlag, Berlin, pp. 137 -152 1981
17 Patterson, c., EI-Sebai, N. A S., A regular boundary method using non-conforming
elements for potential problems in three dimensions, In: Boundary Element Methods in
Engineering. Ed. Brebbia, C. A, Springer-Verlag, Berlin, pp. 112-126, 1982
18 E1-Sebai, N. A S., An Investigation of the Regular Boundary Element Method in Three
Dimensions. Ph.D. Thesis, University of Sheffield 1982
19 Danson, D. 1., BEASY, A boundary element analysis system. In: Boundary Element
Methods in Engineering. Ed. Brebbia, C. A, Springer-Verlag, Berlin, pp. 557 - 557, 1982
20 Aypaci, V. S., Conduction Heat Transfer. Addison-Welsey, USA 1966
Chapter 7
Applications in Geomechanics
7.1 Introduction
The boundary element method is already well established as an efficient numerical
technique to solve a variety of continuum mechanics problems [1, 2]. Solutions
exist for problems as complex as plasticity [3], viscoelasticity [4], viscoplasticity [5]
and other complex time-dependent and non-linear problems. The literature is,
nowadays, very extensive and the interested reader is referred to [2].
In this chapter the boundary element method is applied to solve geomechanical
problems, modelling rock and soil material behaviour. The first model presented
here uses the no-tension criterion and is based on the finite element work of
reference [6]. The criterion was first used by Venturini and Brebbia in the context
of boundary element methods and its numerical implementation is presented in
detail in reference [7]. It consists of considering that the rock medium is unable to
sustain any tensile stresses. The final solution is obtained using an iterative
procedure in which stresses are corrected by introducing a series of initial stress
fields.
Plastic and viscoplastic models are very frequently used to represent the
nonlinear behaviour of soil and rock materials. They have been extensively used
with finite elements since the known classical analysis by Reyes of a circular tunnel
[8]. More recently finite element models have been developed using the more
general Perzina's model [9] to represent viscoplasticity and these models have been
extended to provide plastic and viscoplastic solutions in geomechanics [10, 11]. The
model is applied in this chapter in conjunction with boundary elements following
the original work carried out by Telles and Brebbia [5]. The solution is extended
here by using the overlay technique [12] which allows for a better representation of
material response, as shown in reference [13].
Another important problem in geomechanics is how to model the cracks or
fissures discontinuities often found in the rock medium. As in finite elements [14],
special relationships are used here to model the contact and slip conditions within
a boundary element fomlUlation.
Layers of material, or seams, are also frequently found in soils or rocks. As their
thickness is small by comparison to other dimensions, it becomes inefficient to
consider them in the same manner as other boundary element regions. A special
formulation is described in this chapter which is easy to apply with boundary
elements.
APPLICATIONS IN GEOMECHANICS 143
In the above expression, the two domain integrals correspond to the body force
and the initial stress contributions. The initial stress term, O"~k' is necessary in order
to consider effects such as temperature, shrinkage, swelling and others. The initial
stress integral can also be used to model nonlinear effects such as plasticity, visco-
plasticity, etc. For these cases, an iterative process needs to be carried out. The
fundamental solution used in all cases is the linear Kelvin type solution [15] pre-
sented in the Appendix.
Equation (7.1) can be specialized for an interior point as follows,
Ui(S) = - JUk (Q) Pfk (s, Q) dr (Q) + JPk (Q) Ufk (S, Q) dr (Q)
r r
+ Jbk(q) Ufk(S, q) dQ(q) + JO"~dq) e{mk(S, q) dQ(q) (7.2)
Q Q
After differentiating this equation and applying Hooke's law the stress at an
interior point can be written as [3], [14],
O"ij(S) = - JSijk (s, Q) Uk (Q) dr (Q) + JDijk (S, Q) Pk (Q) dr (Q)
r r
+ JD ijk (S, q) bk (q) dQ (q) + JEijmk (S, q) O"~k (q) dQ (q)
Q Q
1 0 0
8(1- v) [2 O"ij(S) +(1- 4v) O",,(s) bij] (7.3)
Although equation (7.3) only relates to internal points, we will assume that the
boundary stresses have also been included in equation (7.5) since they are usually
required in the solution of nonlinear material problems.
Equations (7.4) and (7.5) are valid for homogeneous bodies. Bodies with region
with different material properties are important in practice. One can analyse them
by dividing the domain into homogeneous subregions.
Consider a domain such as in Fig. 7.1 formed by several different subregions
(Q], Q2." Qi, Q j •. . ). We can start by assuming that each subregion Hi" is an
independent domain and apply equation (7.4),
(7.6)
where the summation notation is not implied.
For the nodes along the interface rij - common to subregions Hi" and ''j'' - the
displacement compatibility and the traction equilibrium are given by,
ui-lJi = 0 (7.7)
pi+ pj= 0 (7.8)
where the superscript represents the subregion to which the displacements or
tractions are related.
Equation (7.6) can also be written for the subregion ''j'' and then imposing the
conditions (7.7) and (7.8) the final system of equations can be obtained by
assembling. This system has the same form as presented in equation (7.4) although
now SOllle of the unknowns displacements and tractions are on the interface rather
than on the external boundary.
After solving the system of equations all external and interface tractions and
displacements are known. Then the stresses at any point can be computed by
applying equation (7.5) separately for each subregion.
The combination of subregions described above presents a particular problem
for comer nodes, where the values of tractions are not uniquely defined.
Traction Discontinuity. In order to understand this problem, consider the two
adjacent elements shown in Figure 7.2. One element is defined by the extreme
points P and R and the other by Sand T where the points Rand S have the same
coordinates.
APPLICATIONS IN GEOMECHANICS 145
(T)
Uz
(PI
p~
"1~":"
Fig. 7.2. Double nodes
X1
Assuming that the stress tensor ~ill be uniquely defined at the double nodes R
and S and writing the traction vector for the adjacent elements as a function of the
stresses, the following condition is derived for the plane case due to the symmetry
of the stress tensor,
(7,9)
where rk) is the direction cosine of element "j" with respect to Xi' Notice that the
Pi components have been assumed in the direction of the global system of
coordinates.
Another condition for the plane case is derived by enforcing the in variance of
the strain trace tensor. Using a linear interpolation function to approximate the
displacements over the two adjacent elements, the invariance condition can be
written,
[(u\P) - U\R» 'l~l) + (U~R) - u~P» 'lP)jlll
+ [( u\T) _u\S» 'l~2) - u~T) u~S» 'IF)j I12
The above two conditions can be introduced into the system of equations
modifying the original equations for one of the corner nodes, This procedure
removes the singularity that otherwise will appear in matrix H and allows different
subregions to be properly combined.
Thin Subregion, The subregion technique presented in this section can solve
systems formed of several homogeneous parts with different material properties,
The technique can also be extended to solve the case of thin layers or seams which
146 APPLICATIONS IN GEOMECHANICS
are frequently present in rock masses. However, in these cases, the interfaces
defined by rock and the seam materials have to be discretized into elements small
enough in order to ensure that nodes on opposite interfaces are at a sufficient
distance from each other in comparison with the size of the elements, otherwise
numerical inaccuracies will arise. As a result the final system of equations obtained
is inconveniently large.
An alternative scheme to improve the boundary element formulation consists of
finding relationships between points on different faces of the thin layer. Consider
the thin layer as shown in Fig. 7.3 and apply one-dimensional stress-strain relation-
ships for shear and compression in a finite difference manner. If the thick-
ness of the seam element, h, is small by comparison with the element length, one
can write,
oW = (u\S) - U\R» Glh
(7.11)
oW (ui
= S) - uiR » EI h
where G and E are the shear and Hooke's moduli respectively and the displace-
ments are taken following the local system of coordinates as indicated in figure 7.3.
Using the traction stress relation one can determine the following equations
relating traction and displacements,
p\R) = (U\R) _ u\S» Glh ,
(7.12)
piR) = (U~R) - u~S» Elh
Two other relationships are provided by the equilibrium conditions i.e.,
p\R) + pIS) = 0 ; piR ) + p~S) = 0 (7.13)
The introduction of these equations into the global system allows for the
representation of a seam.
In all cases the final system of equations is obtained by reordering the matrices
in the usual boundary element manner to take into consideration the boundary
conditions. This gives,
AX= F+ EO'° (7.14)
_ /v
X"u"p,
The vector F takes into consideration the applied body forces and the known
displacements and tractions.
A in equation (7.14) is a non-symmetric matrix which may be banded or fully
populated depending on which way the subregions have been considered. This
matrix is usually well conditioned with dominant diagonal terms.
The solution of equation (7.14) is obtained by inverting A which gives,
X = M + R 0-° (7.15)
where
R=A- 1 E, M=A- 1 F (7.16)
For the stresses (equation (7.5)) we find,
0- = F' + A' X + E' 0-° (7.17)
The above relationships gives the total stresses as a function of the initial stress
field 0-°. If we prefer to write it in terms of elastic stress, equation (7.17) becomes,
o-e = F' - A' X + E* 0-° (7.18)
where
E* = E' - I, o-e = 0- - (l0 (7.19)
Substituting (7.15) into (7.18) gives,
o-e = N + S 0-° (7.20)
in which
s= E* - A'R,
(7.21)
N= F'-A'M
It should be noticed that the vectors M and N represent the solution of a
problem without internal stresses whose effects are modelled by matrices Rand S.
j - 1 --+·1-' 2
9-0-<:>----<)-----<:>--------0----1
-r-- 4-j
~ W
¢-o-o--o----O-------o-----I' f.--o-H-o
1
o Boundary nodes
I
+------
t
- - - - - - - - t ---
- - - Tension
3~, - - Compression
• E,IE,=100 BEM
o E,IE,=l
\
\
\
d_
(f
/
/. Linearly elastic -<Unearly elastic
/ curve // curve
/ //
/
/
.
/
d LL~ng curve
/ /
iea~~ve -------,+=. _ _ . ..
-------+--7 - «tension) E (compression) ~ e(tension)
Unloading curve ~ Unloading curve
Fig. 7.7. Uniaxial stress-strain curve. Path Fig. 7.8. Uniaxial stress-strain curve. Path
independent criterion dependent criterion
150 APPLICATIONS IN GEOMECHANICS
by the equivalent elastic stresses (Je (dashed line in Fig. 7.7). For the path
dependent criterion, the principal stresses are computed using the actual value of
the stresses, (J, obtained after computing the corresponding initial stress field.
Bearing in mind the difference between the two criteria for the verification of
the no-tension condition, the iterative process to model the final solution is illus-
trated by the following steps,
(i) Compute the elastic stress increment for the iteration n. The increment of the
elastic stress is given by (equation 7.20),
(7.22)
in which (n-I)(J0 is the tensile stress obtained in the previous iteration.
(ii) Find the stress vector (Je and (J. They are computed by the following
expressions
(n)(Je = (n-I)(Je + (n)LI (Je
(7.23)
(n)(J = (n-I)(JI + (n) LI (Je
where (n-I)(Jt represents the true stresses obtained in the last iteration.
(iii) Compute the principal stresses. Using (n)(Je or (n)(J as dictated by path
independent or path dependent material behaviour.
(iv) Determine the initial stress vector (n)(J0. The initial stress field to be applied
to the system is computed by taking the values of (n)(J in the tensile principal
direction obtained in (iii).
(v) Compute the true stress vector. The actual stress distribution is evaluated as
follows,
(7.24)
(vi) Verify the convergence. At this stage the convergence criterion must be ap-
plied. If the initial stress vectors are small with respect to the specified
tolerance, the iterative process should stop. Otherwise all steps must be
repeated.
For the two cases studied, the displacements can be computed after applying all
loads or, if necessary at any other stage, using the accumulated initial stress values
in equation (7.15).
The behaviour of rocks for the no-tension criterion described above has been
assumed to be modelled within the context of continuum mechanics. However, it is
know that many types of rocks are continuously disrupted by planar cracks or
fissures. When the rock material is under compressive state of stress, the cracks can
be shut restoring the mechanical continuity of the material, while when under
tension state of stress normal to the fissures, the separation increases and quickly
destroys the existing cohesion. In this case, the rock material loses its continuity
and must be analysed taking into account the possible separation or slip which
may occur between the two surfaces defined by the discontinuity.
In order to model this behaviour using boundary element formulation, the zones
of the continuous material are assumed as subregions, while specific relations
between the values of two opposite nodes along the discontinuity must be defined.
The discontinuity can be interpreted as being formed by special "joint element"
APPLICATIONS IN GEOMECHANICS 151
(Fig. 7.9). The joint element width, h, represents the real gap of the discontinuity
which, in general, can be filled with soft material. The width is taken equal to zero
when the subregions are in contact and the coordinates of two opposite points are
the same. The discontinuity has essentially no, or limited, resistance to a tension
force applied in the normal direction. The slip between the two surfaces of the dis-
continuity is usually assumed to be governed by a linear Mohr envelope which
may be also defined with limited tensile strength.
The equations for the joint can now be written. These equations which are
always dependent on the physical characteristics of the discontinuity must be also
introduced into the final system of equations. They are necessary to model any dis-
continuity problem and are divided in three main groups: (i) the first group
represents the separation condition. This particular condition must be enforced to
the problem when either a real gap between adjacent subregions exists or tensile
stresses in the direction normal to the joint element exceed the admissible limit.
The conditions to be enforced at each node when no soft material is filling the gap
is that the surface traction is zero. (ii) The second group of equations necessary to
model the discontinuity problem is formed by the contact equations,
U\R) - u\S) = 0
(7.25)
p\R) + piS) = 0
i.e. perfect continuity exists. These conditions are valid provided the surface
stresses do not violate a failure criterion. (iii) The last condition to be formulated
refers to the slip between the two surfaces of the discontinuity. Depending on the
stress level along the discontinuity, the two surfaces can slide on each other
without any relative displacement in the direction orthogonal to the discontinuity.
In order to introduce such a behaviour in the system of equations, the following
conditions of displacements and tractions in the direction orthogonal to the dis-
continuity can be written,
(7.26)
U~R) - u~S) = 0
In the direction parallel to the discontinuity, equations relating the shear and
normal stresses can be written according to the assumed Coulomb's envelope. Thus
152 APPLICATIONS IN GEOMECHANICS
100
1
~ 80
E
:;::
.!:::
'C 60 - - theoretical solution
15'
o
•
coarse mesh
fine mesh Coarse mesh Y
I
Fine mesh
40
4 10
R in m
Fig. 7.10. Lined tunnel. Discretization and elastic results
agreement with the theoretical solution. The continuity and the accuracy of the
stress profile on the two discretized interfaces shows the validity of the subregion
technique employed.
The no-tension solution was first achieved by assuming the unlined case in
which the load is applied directly to the rock material. Using 8 linear elements to
discretize both the rock surface and the interface between intact and fissured
material (Fig. 7.11), the final no-tension results are obtained.
The no-tension analysis of the lined case has been carried out using two different
internal cell discretizations and a single boundary mesh consisting of 8 linear
E
" - - Analytical
• Numeric
o Nodes
- - - Cells
8 10
r in m
Fig. 7.11. Unlined case. Discretizations and no-tension results
154 APPLICATIONS IN GEOMECHANICS
140 r - - - - - - - - - - - - - - - - ,
- - theoretical solution
o coarse mesh
120 •
fine mesh
100
-dR(concrete)
~ 80
.!:
/'
/'
~ 60 I
C
Coarse mesh Y Fine mesh
40
2 4 6 8 10
R in m
Fig. 7.12. Lined case. Discretizations and no-tension results
elements and 9 nodes on the internal lining surface and on the interfaces (Fig. 7.12).
The final results for the coarsest mesh seem to be satisfactory; when they are
compared with the theoretical solution the maximum error observed was 4.5%. For
the fine mesh, better agreement between theoretical and numerical solutions was
obtained; the maximum computed error is only 2%.
I
~1====:8=~=2
='8=:1'=4~-j
JI = 0.3
<p=30'
c=O
Plane strain
Fig. 7.13. Circular tunnel and discontinuity surface. Geometry and discretization
APPLICATIONS IN GEOMECHANICS 155
Discontinuity case
Elastic case
\
\
\-----
\
\
\
I
rt
L rt ----------~
t-I~-- 2.64 - - . . ,
Fig. 7.14. Boundary displacements
displacement between the two discontinuity surfaces exists before the application
of the load. Only horizontal loads are applied to the system in order to relieve the
uniaxial system of residual compressive horizontal stresses (all = - 1.0). The
results obtained with the boundary technique are shown in Figures 7.14 and 7.15. In
the first figure the new position of the boundary nodes is presented in comparison
with the corresponding displacements for non fissured case. The second diagram
L
i~ /
r----
r---
V ,- / -
L /
,- ,/
j/
~~'<,.With
A ·th t d·Iscanr·t
!nUl y
/ WI au
/ /
~,-
,- /
./
../
-1
0' 30' 60' 90'
rp-
Fig. 7.15. Normal stress on the tunnel surface
156 APPLICATIONS IN GEOMECHANICS
shows the distribution of the normal stress in the tangential direction over the
tunnel boundary. These results also illustrate the effects of the discontinuity when
compared with the solution obtained using intact rock medium.
7.4 Viscoplasticity
For the analysis of stresses and displacements of any soil or rock structural system
the material properties must be considered as realistically as possible in order to
provide safe and stable design. The previous no-tension and plasticity theories are
based on neglecting the time-dependency of the material properties. However, it is
known that in many geomechanical problems the actual behaviour of the material
is governed by rheological properties. For instance, plastic solutions are acceptable
when the onset of the plastic deformations occurs much faster than the loading
time, i.e. the concept of instantaneous development of permanent and irreversible
strain is valid. On the other hand, the viscous effects become important when
permanent deformations take a comparatively long time to be developed.
In order to introduce time-dependent effects into the boundary formulation, the
Perzyna's model [9] is adopted here. Using this model the stress conditions which
indicate the onset of viscoplastic deformations are governed by a well-known
plastic yield function, i.e. irreversible viscous deformations take place when
(7.29)
where Fis the yield function which can also be represented by,
In equation (7.30) Y (k) stands for the yield stress, k is a history dependent
hardening (or softening) parameter, and !«(1ij) can be interpreted as an equivalent
uniaxial stress value (1e'
In the usual manner for nonlinear problems, the total strain can be separated
into elastic and inelastic parts. In this case in particular, the total strain rate, aij,
can be expressed as,
in which eij and aU stand for the elastic and viscoplastic strain rates.
Adopting the concept of plastic potential from the plasticity theory the
viscoplastic strain rate in equation (7.31) can be written as,
og «(1 ..)
eif= y<([J (FIFo) IJ (7.32)
O(1ij
where y is the viscosity parameter of the material and Fo denotes any convenient
reference value of F, usually Y (k), for the dimensionless representation of ([J.
The function g «(1ij) is also a scalar function of the stress values and equal to
!«(1ij) for associative viscoplasticity cases.
APPLICATIONS IN GEOMECHANICS 157
The function rp is defined from experimental tests and may assume different
forms, as suggested in reference [9] in which the foIlowing representations have
been proposed
rp (FIFo) = FIFo
rp (FIFo) = (FI Fot (7.34)
rp (FIFo) = exp (FIFo) - 1
which correspond to linear, power and exponential types respectively.
After multiplying both sides of equation (7.32) by aij gives,
= Y <rp (FIFo)
.p ) cg (aij)
aij Gij aij ca .. (7.35)
']
Using the general form of equation (7.31) for a continuum problem, the total
strain rate tensor of a particular overlay can be represented by,
(m)"L/].. = (m)i;e.I] + (m)i;p'] (7.40)
where the total strain rate tensor must be the same in all overlays.
Multiplying equation (7.40) by the elastic compliances gives,
(m) . . = (m) ·e._ (m)·p (7.41 )
a" a ,] a'j
158 APPLICATIONS IN GEOMECHANICS
L1 a ij = L1 au - L1 afj . (7.43)
As the same strain pattern is enforced in all overlays the viscoplastic stress
increment tensor, L1afi' is also computed following equation (7.37), i.e.,
(iii) Applying the weighted viscoplastic stress increments, L1 uP, as an initial stress
field, gives the total elastic stress increment vector,
(7.47)
or the corresponding elastic value for each overlay,
(7.48)
(iv) The total stress vector for one overlay "m" can now be obtained at time
t n+ I = t n + L1 t by accumulating the elastic and viscoplastic stress increment
components, i.e.,
(7.49)
then the total stress vector for the whole body is computed once more using
the weighted parameter,
u n+ 1 = (m)u n + 1 h m (m = 1,2, ... , M) (7.50)
APPLICATIONS IN GEOMECHANICS 159
(v) Before computing the next time step, the new value of the equivalent plastic
strain vector can be determined as follows,
(7.51)
in which L18~ is computed employing the equivalent viscoplastic strain rate
given in equation (7.37).
It is important to notice that the steps shown above only model the time
behaviour starting from any state of stress of the body under consideration. In
order to consider the application of any load, whether at the beginning or during
the process, its elastic effects must be instantaneously added to both elastic and
total stress vectors.
The success of the incremental technique shown above for solving elasto/visco-
plastic problems is directly dependent upon the appropriate time step selection.
Then, in order to regulate the time step length, two procedures borrowed from
finite element formulation are adopted here.
The first corresponds to the application of a limit time step formulated for creep
problems [20], and frequently adopted in geomechanics. The time step length is
chosen in order to achieve accuracy in the results and is calculated as a function of
the relation between the equivalent viscoplastic strain rate and the accumulated
equivalent strain.
The technique recommended above cannot avoid instability problems due to
error propagations. Therefore, the time step limit proposed by Cormeau [21] must
also be followed. In this case, the time step length is evaluated based on the
material properties, the adopted yield function and the viscoplastic constitutive
law.
In the examples which follow this section, the viscoplastic algorithm is not only
employed to model time-dependent behaviour but can also be used for plasticity.
As has been known for some time, the plastic behaviour can be obtained by
applying the load in increments and allowing the solution to progress until
stationary conditions are attained. After applying an increment of load the steady
state solution which corresponds to the condition F;§; 0 is approached asymptoti-
cally with the viscoplastic strain rate diminishing as the yield surface is ap-
proached. As for plastic solution, the condition F;§; 0 must be attained during the
whole loading process; the accuracy of such an approach is directly dependent on
the increment lengths adopted.
The final plastic solution was obtained by using the viscoplastic algorithm. The
load was applied in small increments and the convergence was assumed when the
stationary conditions were achieved. The time step and the viscosity parameter
have been chosen in order to obey the stability conditions given in reference [21].
a) Weightless Soil Case. We shall start the analysis of the strip footing by
neglecting the soil weight. The boundary and internal discretizations employed to
run the problem are shown in Fig. 7.16 together with the boundary conditions
prescribed in this case. In common with finite element analysis for this example
(see reference [10]) a closed domain has been chosen. This kind of restriction for
which the displacements are enforced to be zero at points not very distant from the
applied load is a normal procedure in finite element techniques and has proved to
be acceptable for the bearing capacity determinations.
The collapse load predicted by Prandtl was simulated with 1% of error. The load
increment for which the convergence was verified gives the total amount of
21600Ib/ft 2, while the theoretical solution is 21370 Ib/ft2. The problem has also
been solved with the associated Drucker-Prager yield criterion and the same
bearing capacity was obtained.
The load displacement curve for both criteria are shown in Fig. 7.17 where one
can see that the displacements for the Drucker-Prager criterion are larger than the
corresponding values computed using the Mohr-Coulomb surface.
The growth of the plastic zone defined by the Mohr-Coulomb criterion is shown
in Figure 7.18. It began at points underneath the load and only spread all over
the domain when almost all increments of load had been applied.
The determination of the limit load has also been carried out by analysing the
elasto/viscoplastic response. The load is applied in one increment and the final
12ft
25
.10 3
Ib/lt 2 _Pran~'~ulution,,-
'0
20
15
Cl.
"0
a
0
~
10
o
Settlement 01 node A
Fig. 7.17. Load displacement curve
~ p"8.5
OJ] p"11.5
~ p"20
~p"21.6
Load p in 10 3 lblft2
solution is obtained when the stationary conditions are achieved. By solving the
problem several times it was found that a load equal to 21600 lb/ft2 was the largest
value or which the steady solution was obtained. For any load greater than this
value, a continuous growth of the displacements with respect to time was observed.
In spite of being a non incremental procedure, the displacements computed with
this procedure are similar to those obtained following the incremental scheme.
Significant differences have been verified only near the limit load, as shown in
Figure 7.17.
b) Bearing Capacity Considering the Weight of the Soil. This example represents a
more realistic situation if compared with the previous case. For instance, it can
162 APPLICATIONS IN GEOMECHANICS
"'-
"0
C
0
-'
20
15
'w
I
Ii.
10
o 234 5
Settlement of point A
Fig. 7.19. Load displacement curve considering the weight of soil
~ p=10
[[] P =17
~ p=23
~ p=24.2
Load p
in 10 3 lb/ft Z
Unit weight
120 lb Ift3
represent any foundation resting on the ground surface. The self-weight of the soil
y= 120 Ib/ft3 is taken into consideration by assuming a pre-existing stress field in
the material. The vertical stresses are assumed to be equal to the soil weight while
the horizontal values are computed considering the earth pressure coefficient at-
rest, ko, which is given in this case by v/(l- v).
The same analysis carried out for the weightless case has been repeated here. A
limit load equal to 24200 Ib/ft2 has been achieved and is much smaller than the
collapse load computed using the Prandtl's mechanism, which is 25080Ib/ft2•
However, the collapse load assuming Prandtl's mechanism is rigorously an upper
bound. In this case, the ultimate capacity achieved by adopting the Hill
mechanism [22] and assuming valid the Terzaghi superposition [23] is 23450 Ib/ft2,
APPLICATIONS IN GEOMECHANICS 163
only 3% smaller than the numerical result and can be considered a more realistic
solution. Spencer [24] has also achieved a collapse load for this problem by using a
perturbation technique, and his solution, 23800 Ib/ft2, is less than 2% smaller than
the numerical limit obtained. Assuming the Drucker-Prager criterion as for case
"a", the limit load obtained was virtually the same.
The load displacement curve for Mohr-Coulomb yield criterion is shown in
Fig. 7.19 where the theoretical limit loads mentioned above are also presented.
The plastic zone developed during the loading is presented in Figure 7.20. By com-
parison with the first case, the final yield zone is smaller, which can be justified by
the compressive state of stress assumed to exist before the loading.
Parameter Material
Rock Concrete
Before the lining insertion, the opening has an internal diameter equal to 11 ft.
After the application of the concrete support the internal diameter is reduced to
10 ft. The boundary and internal discretizations, together with the geometry of the
problem, are presented in Figure 7.21. As can be seen, only a quarter of the domain
needed to be discretized due to the symmetry. Several loadings are here assumed
to represent different stages of the tunnel construction and they will be illustrated
separately in this section.
a) Unlined Case. This is a situation which often ocurs in practical tunnelling
problems, either due to the deliberate omission of the structural lining or due to a
prolonged delay before the insertion of the support.
Depending on both the viscous properties of the rock and the speed of the
excavation, elastoplastic or visco/elastoplastic solutions are acceptable to model
the stress and displacement distribution around the cavity. Elastoplastic behaviour
is assumed when the time-dependent displacements occur very quickly in com-
164 APPLICATIONS IN GEOMECHANICS
parison with the face advancing. However, the case showing long-time displace-
ments in which inelastic deformations are noticed far behind the face and for quite
a long time is the most common situation.
In order to solve the tunnel excavation assuming visco/elastoplastic behaviour,
Perzyna's model was employed together with an associative Mohr-Coulomb flow
rule. The linear function r/J (see equation 7.34) was applied. The time was defined
in units of lIy (y is the viscosity parameter) and the incremental stepping
procedure has been taken according to the Cormeau's time step. The cohesion and
the friction angle were taken to be 1.44 x 1041b/ft2 and 30 0 respectively. At this
stage, the only load applied is due to the removal of the residual stress on the
opening surface.
Assuming that the load is applied in one single increment, the final elasto/visco-
plastic solution is obtained and the displacements over the tunnel surface are
shown in Fig. 7.22, in which displacements after time equivalent at 5 L1 t are also
presented.
b) Lined Case. In this case, the insertion of the lining is considered in order to
restrain the viscous deformations. The loads, equivalent to the removal of the core,
are applied in one increment and then the lining is placed at a practical distance
from the face. As time and viscous parameter have been hypothetically assumed to
solve this example, the distance from the application of the lining to the face is
considered to be equivalent to a time interval equal to 5 L1 t. Therefore, part of the
viscous displacements have already occurred before the application of the support
(see Fig. 7.22).
The lining is assumed to be of concrete material and will be displaced from the
initial position only due to the viscous deformation of the rock. Figure 7.23 shows
the final steady displacements for the lining, which compare well with the original
finite element solution given in reference [11].
APPLICATIONS IN GEOMECHANICS 165
I
I
I
I
-0- -0- Displ. at 5l\t
I
~ -0-0- Final displ. i -0-0- Lining displ.
1
t ~------
c) Alternative Case. The transference of the stress from the rock to the lining
support obtained in the last example is only due to the time-dependent deforma-
tions which took place in a small region (viscoplastic zone) in the vicinity of the
interface between concrete and rock materials. However, in many tunnel designs it
is convenient to assume that viscoelastic or creep effects occur for any change in
the rock stress level. Using the overlay concept we can define a model to give both
viscoplastic and viscoelastic responses. Here, the viscoplastic effects are modelled
assuming the same conditions adopted in cases "a" and "b", and the viscoelastic
deformations are governed by an equivalent Kelvin-Voight unit. Figure 7.24 shows
the rheological model obtained by representing the rock material using two overlay
models. Also in Fig. 7.24 the parameters adopted to simulate the viscoelastic and
viscoplastic behaviours are presented. The internal discretization for this example
has to be sufficiently increased in order to take into account all significant viscous
deformations.
Overlay 1 Overlay 2
I, =0.5 11 =0.5
E, =7.2·10 7 IbItt 1 E1 =7.2·10 7 Ib/ft 1
P, =0.2 P1 =0.2
[, =0 [1 =2.88·10" Ib/lt 1
1[1, 7], 7]1
1[1,= 0 1[11= 30
r =0.378·1O- 6/month r =0.01 Imonth
d I =0.05868 month
'::Tc.:::Q~~
;1'=··. ~
7.5 Conclusions
In this chapter the formulation of the boundary element method to deal with the
no-tension and discontinuity processes of geomechanics has been presented. The
examples presented here show the accuracy of the solutions. Boundary elements is
attractive for solving geomechanical problems because of the considerable reduc-
tions in the data required to run a problem and the possibility of modelling
properly boundaries extending to infinity.
The chapter explains how to formulate the nonlinear algorithm required for the
solution of no-tension materials. They have also been extended to deal with
viscoplastic problems using the overlay technique and Perzyna's model. The
procedure can be used as well for the solution of time independent plasticity as a
particular case.
In the section 2 of this chapter the boundary integral equations for displacement
and stress determinations have been introduced. They involved the, so called,
Kelvin's fundamental solution and other tensor forms. All these tensors correspond
to responses at a source point "s" due to a unit load applied at a field point "q" in
an infinite plane domain.
For completeness we resume here all of these forms,
fundamental solution for displacements:
utes, q) = 8 n (1-
-I
v) G
{(3 - 4v) In r bij - r,i r,j} (AI)
P0, (s, q) = 4 n (1 I- or
{ [(1- 2 v) bij + 2 r,i r,jl --;- - (I - 2 v) (r,i Y/j - r,j Y/i) } (A2)
v) r un
etmk (s, q) = 8 n (I ~ v) G r {(l- 2) (r,k bim + r,m bik ) - r,i bmk + 2 r,i r,m r,k] (A3)
2G
Sijk (s, q) = 4 n (1- v) r2 {2 r,n [(1- 2 v) bij r,k + v (bik r,j + bjk r,;) - 4 r,J,j r,kl
+ 2 v (Y/i r.j r,k + 'lj r,i r,k) + (1- 2 v) (2 Y/k r,i r,j +
+ 'lj bik + Y/i bjk ) - (1- 4 v) Y/k bij } (A4)
I
D'k (s q) =
lJ' 4 n (1- v) r
{(I- 2v) [bik r) + b'j'k r i- bi/,r kl
","
+ 25 ir),r
, ,
k} (A5)
168 APPLICATIONS IN GEOMECHANICS
1
Eijmk (s, q) = ? {(l- 2 v) [(lik (ljm + (ljk (lim - (lij (lmk + 2 (lij r,m r,k]
4n(l-v)
+ 2 V [(lim r,j r,k + (ljk r,i r,m + (lik r,j r,m + (ljm r,i r,k]
+ 2(lmk r,ir,j- 8r,i r,jr,m r,k} (A. 6)
In all the above expressions the r,i represents the derivatives of the distance be-
tween "s" and "a" with respect to the Xi Cartesian coordinate of "q", i.e.,
r·=
or (s, q) = Xi (q) - Xi (s)
,I OXi (q) r (s, q)
where
References
1 Brebbia, C A, The Boundary Element Methodfor Engineers. Pentech Press, London 1978
2 Brebbia, C A, Telles, 1. C F., Wrobel, L. C, Boundary Elements Techniques. Theory and
Applications in Engineering. Springer-Verlag, 1984
3 Telles, 1. C F., Brebbia, C A, Plasticity. In: Progress in Boundary Element Methods. Ed.
by Brebbia, C A, Chapter 5, 2, Springer-Verlag, 1981
4 Rizzo, F. 1., Shippy, D. 1., An application of the correspondence principles of linear
viscoelastic theory. SIAM, 1. Appl. Math., 21, pp. 321- 330, 1971
5 Telles, 1. C F., Brebbia, C A, Elastic/viscoplastic problems using boundary elements.
Int. 1. Mech. Sci., Vol. 24, pp. 605-618,1982
6 Zienkiewicz, O. C, Villiappan, S., King. 1. P., Stress analysis of rock as a no-tension
material. Geotechnique, 18, No. I, pp. 56-66, 1968
7 Venturini, W. S., Brebbia, C A, Boundary element formulation to solve no-tension
problems in geomechanics. In: Numerical Methods in Geomechanics. Ed. by Martins, 1. B.,
NATO ASI Series, REIDEL, 1982
8 Reyes, S. F., Deere, D. U., Elastic-plastic analysis of underground openings by the finite
element method. Proc. 1st Congress of the Int. Soc. of Rock Mechanics, 2, pp. 477 -483,
Lisbon 1966
9 Perzyna, P., Fundamental problems in viscoplasticity. Advances in Applied Mechanics, 9,
pp. 243-377, Academic Press, New York 1966
10 Zienkiewicz, O. C, Humpheson, C., Lewis, R. W., Associated and non-associated visco-
plasticity and plasticity in soil mechanics. Geotechnique, 25,671 -689. 1975
II Zienkiewicz, O. C, Humpheson, C, Viscoplasticity: A generalized model for description
of soil behaviour. In: Numerical Methods in Geotechnical Engineering. Ed. Desai, S. et aI.,
McGraw-Hill 1977
12 Owen, D. R. 1., Prakash, A, Zienkiewicz, O. C, Finite element analysis of non-linear
composite materials by use of overlay systems. Computers and structures, 4, pp. 1251-
1267, 1974
13 Venturini, W. S., Brebbia, C A, The Boundary Element Method for the Solution of No-
tension Materials. In: Boundary Element Methods. Ed. Brebbia, C A., Springer-Verlag,
Berlin Heidelberg New York 1981
14 Goodman, R. E., Taylor, R. L., Brekke, T. L., A model for the mechanics of jointed rock.
1. Soil Mech. Found. Div., ASCE, 94, pp. 637-657, 1968
15 Love, A E. H., Treatise on the Mathematical Theory of Elasticity. Dover 1944
16 Hocking, G., Stress analysis of Underground excavations incorporating slip and separation
along discontinuities. In: Recent Advances in Boundary Element Methods. Ed. Brebbia,
C A, Pentech Press, 1978
17 Jaeger, 1. C, Cook, N. G. W., Fundamentals of Rock Mechanics. John Wiley and Sons,
Inc. New York 1976
APPLICATIONS IN GEOMECHANICS 169
Applications in Mining
8.1 Introduction
Mining is the term used for the extraction of mineral from the earth's crust. The
problems of analysis which arise are thus almost exclusively concerned with the
excavation of ore from the "host material" and the stresses and displacements
induced by the consequent mining excavations. The determination of these stresses
and displacements then allows predictions to be made about the stability of the
mine working and may affect the design of the mine layout. There are basically
two ways in which ore is deposited, namely,
1. Tabular: In seams, lenses or strata of large extent in two dimensions but with
relatively small thickness in the third dimension.
2. Massive: A three dimensional body of complicated shape.
For tabular orebodies the terms "roof' and "floor" or "hanging wall" and
"footwall" are used for the surrounding rock mass surfaces adjacent to the seam.
Figure 8.1 shows a cross-section of the Mount Isa Mine, Queensland, in Australia
with both tabular and massive orebodies. One method of excavation is to create
openings (or stopes) separated by ribs of the ore material to form pillars. These
pillars provide support and prevent roof or floor collapse. The design or
dimensioning of these pillars is critical to the efficiency of the mining operation.
In a second approach used mainly in coal mining, the roof is kept stable only in
the working area. Collapse is permitted in mined areas away from the current
excavation zone. The design of the support systems in the working area is vital to
the success of this type of operation.
For massive orebodies, or steeply inclined tabular deposits, ore can be extracted
by a method called "open stoping" in which large voids (or stopes) are formed by
blasting. The dimensions of these "stopes" can be very large (in the order of
100 m). After blasting, the fragmented ore is drawn from the base of the stope. The
stability of the walls of the stope (hanging wall and footwall) is of concern mainly
because the collapse of host rock with no mineralisation into the stope results in
dilution of the ore causing inefficiencies in the mineral processing operations. In
many cases, cable dowels or rock bolts are installed to attempt to prevent hanging
wall collapse. In the "open stoping", a certain pattern is followed leaving enough
material in the pillars to ensure regional stability (Fig. 8.2). After the ore has been
withdrawn, the stopes are backfilled with material which may be a mixture of
APPLICATIONS IN MINING 171
West East
- - 8 level 328 m
- - 15 level 728 m
Fig.8.1. Cross section of Northem Mine area, Mt. Isa, Queensland, Australia
Orebody continues
to the south
-8-- N
Filled cemented
Producing hydrau lic titt
O Empty or filling
1:71 Filled cemented
~ rockfitt
250 m r.:ol
l!tiJ Filled, rockf ill
D
DeSigned
Fig. 8.2. Plan of open stoping in the 1100 ore body at Mount Isa
172 APPLICATIONS IN MINING
milltailings, cement and possibly rock fragments (rock fill). Extraction of the pillar
ore follows. The objective is to achieve maximum extraction of ore with a
minimum dilution and at the same time maintaining regional stability.
The above examples of mining methods are given to highlight some of the
problems encountered in a mine operation. It can be seen that the analysis
problems associated with a mine are large, complex and three-dimensional (see for
example the large number of stopes and pillars in Figure 8.2). To contain the
analysis effort to reasonable proportions, simplifications have to be introduced and
these have influenced the development and application of Boundary Solution
methods.
The papers by Salamon [1-4] seem to be the earliest references on the application
of Boundary Solution procedures to mining problems. Salamon concentrated on a
solution for the problem when the thickness of the reef or seam is negligible
compared with its overall dimensions. A good approximation to this problem is to
assume that the thickness of the excavation is zero but that the displacements of
the roof are different from those of the floor. A fundamental solution for a
displacement discontinuity or dislocation in an infinite medium is derived and a
"Face Element" approach used to obtain the solution for the actual excavation
problem. The "Face Element" method is, in fact, identical to the Boundary
Element method with constant Boundary Elements. The approach by Salamon was
generalised and adapted to digital computers by Starfield [5] and Crouch [6]. The
method has become increasingly popular in recent years and many workers have
contributed to the further developments which have resulted in improvements in
capability and economy [7 - 10]. The popularity of the method with mining
engineers stems from the fact that the data input is extremely simple. Many
problems which are too complex and too large to analyse by a full 3-D analysis can
be solved.
Deist and co-workers [11-12] were the first to propose a method for analysing
excavations in massive ore bodies. In this method, the excavation surface is
discretised into "Face Elements" and a distribution of singularities is determined
which reduces the traction at the center of each element to zero. The analysis is
fully three-dimensional and the Kelvin solution is used. The method is essentially
an indirect Boundary Element method (with constant Boundary Elements) because
the displacements are not determined directiy but rather from the distribution of
singularities. In the implementation on the computer, Deist has been particularly
concerned with efficient and economic solution of the system of equations.
Complex mine configurations in some instances required discretisations of up to
4800 Face Elements with 14400 degrees of freedom. The solution of such large
systems of equations is a formidable task especially when it is realized that the
coefficient matrix is fully populated and non-symmetric. Deist uses an accelerated
APPLICATIONS IN MINING 173
1) Subscripts denote direction (i = 1,2, 3;) = 1,2,3). For example, Tij is the
traction component in Xi direction due to a source in Xi direction. The usual
summation rules apply.
2) Values in parenthesis denote location. For example, Tij (x, y) is the traction at
location x having the co-ordinates Xi due to a source at location y, co-ordinates
Ji. Alternatively, for constant Boundary Elements, aij(m, m) is the influence
coefficient at element m due to a source at Element n.
174 APPLICATIONS IN MINING
For the mined portion of the seam all components of stress must be reduced to
zero. That is, the induced stresses at the excavation surfaces must exactly balance
the premining stress components P13, P23,P33. Thus,
an =- Pi3 (8.2)
In the unmined portions the following must hold if the seam material is infinitely
stiff:
di=O (8.3)
For an elastic seam the following relationship between stress and displacement
exists
a33 = Eft d33 (8.4)
ai3 = Glt dn (i oF 3) (8.5)
where E and G are the compression and shear moduli of elasticity of the seam and
tis a fictitious thickness.
Solutions for the excavation in seams with any geometry can be found by sub-
dividing the region of interest into a large number of small square zones or
Boundary Elements of area a. The displacement discontinuities and stresses are
APPLICATIONS IN MINING 175
In the above, d;(m) is the average dislocation at element m and the summation is
for all elements M The coefficients aij(n, m) are fundamental solutions for
influences of unit dislocations. The fundamental solutions have been derived using
harmonic functions by Salomon [1] and Starfield [5]. The solutions can also be
obtained by using the Somigliana identity and the Kelvin solution (see Appendix B).
In the early displacement discontinuity programs, a regular matrix of elements
was assumed and an element defined by its row and column number. In this
regular matrix excavated elements were flagged within the grid and any excavation
pattern could be reproduced. The use of a regular grid allows the influence
coefficients to be a function of row and column numbers (or more specifically
differences in row and column numbers) thus making their evaluation more
efficient. In addition, the data input is Simplified to a great extent. For mine size
problems however, the number of unknowns tends to become large and in the
computer programs developed to date an iterative solution method with "lump-
ing" has been used.
The explanation of the method herein has been restricted to the excavation of a
single seam in an infinite domain. The formulation can be extended to excavations
of more than one seam in proximity to each other in either infinite or semi-infinite
media. These extensions are discussed in detail by Sinha [26]. When an iterative
solution of the equations is used it is also possible to enforce Mohr-Coulomb and
no-tension conditions in the unmined portions. Here the yield conditions are
checked every few iterations and the stresses adjusted as required. This process
seems to work [7] although results from non-linear analyses have not been tested
rigorously.
z~
y
Fig. 8.4. Mesh of displacement discontinuity elements. (Courtesy of Mt. Isa Mines Ltd.)
176 APPLICATIONS IN MINING
level:
13E
138
13/l
"'~
-!:; E
(/) "-
=""
~ ::E
Cl.
w
=
31
ZOO m
Figure 8.4 shows typical grid patterns used for the analysis of tabular seams. The
pillar stresses for particular excavation stages at Northern 5 orebody at Mount Isa
[68] are shown in Figure 8.5. The values of pillar stresses are found to be in good
agreement with measurements. While the program gives good results for average
pillar stresses, the detailed stress components in the vicinity of the edges of the
excavation are found to be unreliable. Better results can be obtained by assuming a
linear or parabolic variation of dislocation and stresses instead of constant or
average values. Development in this direction has been carried out by Diering [27]
and more recently by Crawford [10]. Displacement discontinuity programs have
APPLICATIONS IN MINING 177
been extensively used by mining engineers mainly in South Africa, Australia and
the USA. The reason for their popularity is the ease in data input which is
oriented towards the methodology used in mining. Also, the solution of thin
tabular excavations could prove to be difficult with classical Boundary Element
methods because of the proximity of the roof and floor surfaces, a situation which
can lead to ill-conditioning of the coefficient matrix.
S (k) Element k
Fik)
Pi!)
Element I
Fig. 8.6. The direct Boundary Element method
178 APPLICATIONS IN MINING
In equation (8.11), Tij (x, y) is the fundamental solution for the tractions at x due
to a source at y and y(l) is the centre of element t. For t = k i.e. when the source is
at Element k the fundamental solution tends to infinity and the integral does not
exist. For smooth surfaces Deist [12] has shown that the coefficient can be
computed from:
(8.12)
where (jij is the Kronecker Delta.
Equation (8.10) can be solved for the fictitious forces Pj (l). The displacements at
point x inside the continuum can be computed from:
L
Ui(X) = L. A Uij(x, l) Pj(l) (8.13)
where 1=1
200
400
600 ~
a;
Ll
oS
Cl.
'"
Cl
1000
m
1200
Fig. 8.7. Cross section of a mine involving multiple pillars
A direct Boundary Element formulation was first proposed by Rizzo [13] in 1966.
In this formulation Betti's reciprocal theorem is used to obtain an integral equation
on the boundary in terms of displacement (u;) and traction (t;) values. This leads
to:
Cij(X) Uj(x) + JTjj(x,y) Uj(Y) dS= JVij(x,y) tj(Y) dS (S.lS)
S S
where Tij and Vij are fundamental solutions for the tractions and displacements
and cij (x) is t (jij for smooth boundaries. The integral equation may be solved by
approximating the boundary variables Uj, tj as piecewise constant over a finite
area. Thus the surface of the problem to be analysed is divided into M Boundary
Elements and (S.IS) rewritten in algebraic form as:
M M
cijUj(n) + L.
m~1
LJTij(n,m)uj(m) = L.
m~1
LJVij(n,mH(m) (S.19)
where
LJ Tij(n, m) = J Tij(x(n), y) dS(y) (8.20)
and SCm)
where x(n) is the location at the centre of Element nand SCm) is the surface of
element m. In three dimensions the number of equations in (8.19) is 3 M. For an
180 APPLICATIONS IN MINING
II = L LJ Ulj(n, m) tJ(m)
m~1
M
(8.23)
12= L LJU2j (n,m)tJ(m) etc.
where m~1
(8.24)
and rfij is the premining stress tensor.
The coefficient matrix [A] is fully populated and unsymmetric. The displace-
ments at the boundary are determined by solution of (8.22). The displacement at
an interior point x is
M M
Ui(X) = L L1Uij(x,m)tj(m)- L
m~1 m~1
L1Tij(x,m)uj(m) (8.25)
where the third order tensors L1 Dikj and L1 Sikj are determined from differentiation
of (8.25).
The first application of the direct formulation to mining problems was reported
by Cruse [14], who analysed a room and pillar structure representative of the
excavation geometry of a coal mine. A more recent application of the direct B.E.M.
with the displacements and tractions assumed constant over each Boundary
Element has been presented by Diering [35]. Figure 8.8 shows the geometry used to
model the interaction between an open pit and underground excavations. A plane
of symmetry was assumed thus reducing the discretisation effort by one half. For
this problem 350 Boundary Elements were used and it is remarkable that the
solution was achieved on a 16 bit Data General mini-computer. This was made
possible by using a lumping procedure and solving the system of simultaneous
equations iteratively. The results of the analysis high-lighted areas of high tension
or high shear on specified planes.
Ricardella [36, 37] appears to be the first to suggest a higher order variation of
displacements and tractions across the Boundary Elements. He implemented a
linear variation in a two-dimensional Boundary Element program. This program
has been developed further and applied to the two-dimensional analysis of mining
problems by Crotty and Wardle [38]. A 3-D analysis program using triangular
Boundary Elements and linear variation of the unknowns as developed by Cruse
and modified by Wardle to handle excavation problems but has not yet been
applied to large scale mining problems [39]. Lachat and Watson [15, 16, 34,40,41]
have written programs which use curved isoparametric elements with parabolic
APPLICATIONS IN MINING 181
Plane of
symmetry
./
Crown pillar Underground
excavation
Pillar
Stope
For a node where F T > 0 separation occurs and shear can no longer be
supported. The excess tractions tj are
t~ = - ts and t~ = - FT (8.33)
For a node where Fs> 0 inelastic slip occurs at the interface. The excess traction
parallel to the interface is
t~ = - sign (t s ) . Fs (8.34)
The non-linear problem may be solved iteratively by redistributing excess tractions
at each step. The iteration continues until the yield conditions are satisfied at all
points of the interfaces. The procedure has been described for 2-D problems here
but extension to 3-D involves only one additional tangential traction component.
The yield and flow conditions can be written to include peak and residual strengths,
joint dilatancy and other features which give a more realistic description of rock
fractures. The treatment of non-linear material properties is identical to the one
used for Finite Element and there is ample literature on the subject [43].
For rigid interfaces, the coefficient matrices [A] for each subregion are
assembled into a global matrix [A]G as shown in Fig. 8.9 for a problem consisting
of two subregions. It can be seen that the matrix is no longer fully populated. Thus
the effort in establishing and solving the system of equations can be significantly
reduced when subregions are used. The use of subregions, therefore, is also
advisable in cases where the material is homogeneous. Watson has used this
technique extensively for three dimensional stress analysis problems [34]. Special
solvers for the solution of such sparse matrices have been developed [44]. An
application of the multi-region Boundary Element method to a mining problem
was undertaken by Hocking [44]. He analysed the extraction of a massive orebody
at Mount Isa Mines (Queensland, Australia) using a highly idealised model
(Fig. 8.10). Here each pillar (shaded) is treated as a subregion and symmetry con-
ditions have been applied in the x-y, x-z and y-z planes. The mesh consists of
1-----1------1- - - - -
Zero
coefficients
_ . _ _ . _ '---_-1- _ _ _ _ _
~
55m A
/s~H
x
Fig. 8.10. Idealization used for the 1100 orebody at Mount Isa (from [46])
seven regions with over 100 isoparametric Boundary Elements of parabolic shape.
The main result obtained by the study were the pillar stresses at different stages of
excavation. However, these are of little practical use since the actual problem is
much more complicated than the simple idealisation and the boundaries used. One
obvious limitation must have been the number of Boundary Elements needed to
discretise each of the 30 excavations in more detail. In three-dimensional stress
analysis by the Boundary Element method, the costs of constructing and solving the
system of equations escalates rapidly with the number of excavations making the
solution of problems involving the thousands of elements impossible unless special
techniques are used. This is discussed in more detail later.
When the infinite medium is divided into subregions the interfaces theoretically
extend to infinity and an infinite number of Boundary Elements are needed to
discretise the interface. Such problems can be solved by truncating the interface
discretisation at points which are sufficiently far away from the region of interest.
This method has been applied successfully by Crotty [46] and Bolteus [47]. A better
and more accurate method might be to use infinite Boundary Elements as
suggested by Watson [34].
Rigid-plastic interfaces can be analysed by redistributing excess tractions
iteratively. This process may be included in the iterative solution of the system of
equations. An example of an analysis of underground excavations in faulted rock is
shown in Figures 8.11 a and 8.11 b. The figures show results of a program BITEMJ
[46] which has been recently released by CSIRO, Division of Applied Geo-
mechanics.
Another way of analysing excavations in faulted rock is by a combination of
Boundary and Finite Element methods. In this approach two or more Boundary
Element regions are connected with each other by joint Elements. The coupled
analysis procedure will be discussed later.
APPLICATIONS IN MINING 185
700~------------~----------~ ,-------------,------------,
m
I I \ x "'f.. "i-
+ 'f- 'f.-'f.-
-\- 'j.. --....- -
... ---X+"i-'I-
-XD'
600 X --r~---- ----+-/x-+'j.. 'l<
f x ---+----r-/'---..,..'x, '<
0
~-
'i I ",-
\
\ \ 1/1 +
;/ + ;/
"-
"J.....
'" "-
......... -.,.... ____ -+- ------xxx
;<
In addition to major faults, pre-existing fractures and fissures are present in the rock
mass. These are classified under the term "structural discontinuities" and in most
cases, their frequency and continuity is such that it is uneconomic to model them
individually.
One approach has been to consider the combined or overall effect of these
structures i.e. consider the rock mass as an elasto-plastic material. Mohr-Coulomb
or Drucker-Prager yield conditions have been applied to materials which are
essentially isotropic, i.e. where the structural discontinuities are distributed in a
random manner. In cases where fractures have a predominant direction, a multi-
laminate model [48] can be used. In this model, Mohr-Coulomb and no tension
conditions are checked on certain planes only with stress redistribution occuring in
each plane. This model is in many respects similar to the joint model discussed
earlier except that the fracture is not modelled explicitly. Instead, the stresses at a
number of interior points are checked. For a point at which the failure criteria are
exceeded, the excess stresses are treated as initial stresses and redistributed. This
initial stress approach is well known to workers in the Finite Element method [43].
It has been shown recently that these techniques can also be used with the
Boundary Element method but require a discretisation of the continuum in the
zone where elasto-plastic straining occurs. This discretisation is needed because the
redistribution of the excess stresses involves the computation of a volume integral.
Equation (8.18) is rewritten for elastoplastic problem as:
cij (x) Uj(x) + S Tij (x, y) Uj (y) dS = S Uij(x, y) tj (y) dS + SEijk (x, y) aJk (y) dV
s s v
186 APPLICATIONS IN MINING
where Eijk is the fundamental solution for the strain at x [49]. To compute the
volume integral the interior is divided into Finite Elements or cells. In the simplest
case, rectangular cells are used and the stresses assumed to be constant over the
cell. The discretised form of (8.35) is
M M L
CjUi(n) +I ATij(n,m)Uj(m)= I AUij(n,m) tj(m) +I AEijk(n,l)aJk(l)
m=l m=l i=l (8.36)
where L is the number of cells and
AEijk(n, /) = S Eijk(x(n),y) dV(y) (8.37)
V(l)
. Ground surface\
H=124 m y
~~~----.~----
x
Fig. 8.12. Boundary discretisation and internal Fig. 8.13. Elastic and no-tension re-
cells of cavern. (Courtesy of C. A. Brebbia [50]) sults. (Courtesy of C. A. Brebbia [50])
APPLICATIONS IN MINING 187
The process continues until the yield conditions are satisfied at all points. The
elasto-plastic analysis thus involves an interior stress computation for each cell
centroid, a volume integration and a solution for each iteration.
Examples of non-linear analysis in geomechanics are given by Venturini and
Brebbia [50] and in Chap. 7 of this volume. In Fig. 8.12, the discretisation of an
underground rock cavern into Boundary Elements and internal cells is shown. The
result of the analysis assuming the rock to be a no-tension material is shown in
Figure 8.13. There has been an increasing amount of research into the application
in plasticity. However, elasto-plastic Boundary Element analyses have, to the
authors' best knowledge, not been applied to mining problems.
Interface
no des Free
,.---A--, / n odes
!
Boundary
/ elements
\
\\
Finite elements
Fig. 8.14. Coupled finite element-boundary element mesh
displacements and tractions at the interface are {a}n and {b}n and the correspond-
ing values at the free surface are {ah and {bh. The Boundary Element equations
can be written as an extension of equation (8.22).
where {ah and {h}n are the solutions with the interface nodes fixed and the
columns of [Q]I and [Q]n are solutions for unit values of displacements at interface
nodes. The tractions at the interface can now be expressed as a function of the
displacements at the interface by
71:=2I Su;t;dS-
- S-
ujtjdS (8.43)
s s
APPLICATIONS IN MINING 189
where i are tractions due to displacements and i other tractions. For the interface
the coupled equation is obtained from this (see ref. [53]):
([K]F.E. + [K]B.E,) {a}n + {F} = 0 (8.44)
[K]B.E. is the stiffness matrix for the interface nodes. This matrix is symmetric and
for two-dimensional problems has the size of 2Mx 2M where M is the number of
interface nodes. The vector {F} is the nodal force vector. After the stiffness matrix
and the load vector have been obtained for the interface, the analysis proceeds
using standard Finite Element software. The Boundary Element region is treated as
a Super Finite Element and assembled in the usual manner.
When non-linear problems are analysed, the iterations now only involve degrees
of freedom of the Finite Element nodes including the interface. In the non-linear
analysis by the Boundary Element method, the degrees of freedom of all the
boundary nodes have to be considered. Also, in the Finite Element analysis only
the nodes of each Finite Element contribute to the stress computation at points
inside it. This is in contrast to the non-linear solution by the Boundary Element
method where each interior stress computation involves all the Boundary Element
nodes. For problems where the number of interface nodes is small compared with
the total number of boundary nodes a coupled analysis is thus expected to be more
efficient. Some typical non-linear mining problems analysed using the coupled
Finite-Boundary Element program BEFE [69] are shown herein.
Example I. Mine Pillar Problem [24]. In Fig. 8.15 a problem is shown where ex-
cavations have been made in virgin ground. As a result, the material between the
excavations, termed a "mine pillar", is subjected to high stress. This situation
occurs frequently in mining especially in coal mining although in particular
applications the dimensions will vary and several pillars may be present. Because
of the high stress, the pillar is expected to behave non-linearly and thus it has been
discretised into Finite Elements. In the example, the pillar was assumed to have
planes of weakness inclined 30 0 to the horizontal. These planes were assumed to
have an angle of friction of 10 0 and a cohesion of 1 MPa. These properties were
modelled using a viscoplastic multi-laminate model [48]. The vertical stresses in
the pillar for the elastic analysis are shown in Fig. 8.16, together with those for
the first 3 time steps in the visco-plastic analysis. The vertical pre-mining stress
was taken as 30 MPa.
Example II. Crown Pillar Analysis [54]. This is similar to the previous example
except that the excavation is inclined to the horizontal and the shape of the
material between the excavations is more complex (see Figure 8.17). Because of the
Parabolic boundary
elements
/
.' .,
- Parabolic finite f
1"1
" elements 20m
~
I. 90 m 90m
.1
Fig. 8.15. Mine pillar problem: mesh and dimensions
190 APPLICATIONS IN MINING
150.-------.-------.-------.-------~
MN/mz
~ ~.
~100~------~--------~------_,'------~
en
c:
.3
o
V>
V>
~
V>
0L-______J -_ _ _ _ _ _- l_ _ _ _ _ _ _ _L -_ _ _ _ _ _
~
10 o 5 m 10
Distance from centre
Fig. 8.16. Pillar stresses for the first 3 time steps for pillar with planes weakness at 30 0
Fig. 8.18. Crown pillar problem: displaced Fig. 8.19. Crown pillar problem: visco-
shape plastic zones
Finite
elements
the layers of Finite Elements. The angle of friction has been assumed to be 10°
with no cohesion. This value corresponds to particularly smooth and graphitic
bedding planes. Figure 8.21 shows the displacements which occur in the hanging
wall after the excavation has been made in the initially stressed rock mass. It can
be seen that there is significant differential slip occuring near the abutment and,
associated with this, an opening of fractures. This is an example of the coupled
analysis where the problem could not have been analysed as efficiently by either
Finite or Boundary Element methods.
rule, two Boundary Elements should not be closer to each other than half of their
length otherwise the solution can become unstable}.
There are some problems in mining that involve both massive and seam type
excavations and there is obviously some benefit in having both formulations in one
computer program. The two methods are in fact very similar except that instead
of displacements, the closure and rides are the variables in the displacement
discontinuity method. In the Boundary Element equations, the influence of the
displacement discontinuity element has to be included i.e., Eq. (8.19) is rewritten:
M L M
cijUj(n} + ~ LlTij(n,m} Uj(m} +~ LlTij(n,l)~(l)= ~ L1Uij(n, m} tn(m} (8.45)
m=1 1=1 m=1
Present
underground mining 1000m
Section A-A
Fig. 8.22. Schematic representation of open pit and tabular underground excavations (from
[36])
Y·10000E
-50
X.21000 N
The equations (8.45) and (8.46) are then solved for the values of displacements and
displacement discontinuities. For a seam inclined to the global axes, the tractions
and displacements have to be converted to a local coordinate system. Details are
given by Diering [35].
Figure 8.22 shows a typical application of the mixed method ~here a pit
excavation interacts with a tabular excavation. This problem was analysed by
Diering [35] who discretised the Pit surface into Boundary Elements and the seam
into Displacement Discontinuity elements. A typical result is shown in Fig. 8.23
which depicts the contours of displacement of the ground surface. Similar mixed
procedures have been used by Hocking [21] to analyse excavations with fault zones.
Acknowledgements
Some of the work presented in this chapter has been sponsored by Mount Isa
Mines. Ltd., Mount Isa, Australia. The permission to publish mine cross-sections
of the Mount Isa mine and the results of a mining problem analysed by the
displacement discontinuity program NFOLD is gratefully acknowledged. Thanks
must go to T. Diering for permission to include some figures from his thesis here.
The permission of the CSIRO, Division of Applied Geomechanics to include
results from program BITEMJ is acknowledged. Finally, the authors would like to
thank C. A Brebbia for making available plots of the no tension analysis.
196 APPLICATIONS IN MINING
All Boundary Element formulations discussed here lead to large systems of simul-
taneous equations. These coefficient matrices are fully populated and non-sym-
metric. This contrasts with the sparsely populated and symmetric matrices of the
Finite Element method. Thus more storage is required and more time used in the
solution of a given number of equations. On the other hand, however, the number
of unknowns is much less in a Boundary Element discretisation.
Basically, two different approaches can be used for the solution, that is, iterative
and direct methods. In addition, hybrid schemes are possible. One quite popular
iterative solution technique is the Successive Overrelaxation technique.
The solution of
[AJ {a} = {b} (8.47)
TableS.1
Gauss S. O. R.
An inspection of the matrix [A] shows that for most problems many of the off-
diagonal coefficients are very small compared with the diagonal terms. This is
because the Kernel functions decay rapidly with distance, particularly for the
three-dimensional formulation.
One way of exploiting this property of [A] has been proposed recently by Bettess
[57]. In this method, the smaller off-diagonal terms are replaced by zero and a
sparsely populated matrix obtained. This matrix is then solved using a skyline
solver which eliminates most of the operations on zero coefficients. The results are
then multiplied with the full matrix [A] and compared with the original right hand
side. The difference between the results obtained is then taken as a new right hand
side and the system resolved using the factorised sparse coefficient matrix. The
process is a combination of both the direct and iterative method. For the problems
considered in [58], the solver reduced the execution time by a factor of about two.
This factor could increase significantly for very large problems, but no further
details are available.
Another method, perhaps more physically obvious, has been developed by Deist
and used by others. Called the "lumping" method, it takes advantage of the fact
198 APPLICATIONS IN MINING
that inter-element influences decay rapidly with distance. However, these remote
influences can not be neglected because their combined effect is still significant.
For constant Boundary Elements, the lumping procedure is as follows: If an
Element is sufficiently far away from a group of Elements, its effect can be
evaluated on a central point within the group. Similarly, the effects of a group of
elements on any remote element can be represented by the influence of a "com-
posite" element located at a central point of the group. It can be seen that this
method will introduce sparseness into [A] because for distant elements only one
point in the lump contributes to the coefficient matrix. The bigger the lumps the
sparser the matrix will become. In addition, fewer coefficients will have to be
calculated. To make the "lumping" procedure fully automatic is not straight
forward and details of the method are discussed by Deist [12]. Diering [35] has
recently adopted a much simpler approach where lumping information is input by
the user. That is, the user defines "lumps" of elements which are used to evaluate
coefficients for distant elements. In the "lumping" procedure, the matrix [A] will
become sparsely populated but not banded. The iterative solution techniques are
not affected much by the lack of banding, but for the direct solution, efficient
solvers for unbanded sparse matrices are not freely available. For this reason, the
"lumping" has only been used to date in conjunction with iterative solvers.
However, there are no technical reasons why lumping can not be implemented into
direct solvers.
x,
and c( = (1 + v)
4n E (1- v)
(8.57)
C2 = 3 - 4v.
In equations (8.52) and (8.53) the elastic constants are E, Young's Modulus and v
Poisson's ratio. The stress at x is given by
(8.58)
with
(8.59)
In the above, A and G are the Lame constants. Differentiation of equation (8.55)
and substitution in equations (8.59) and (8.58) gives the stress components at x due
to a unit source in directionj at y.
C3
DUk (x,y) = 2R2 [C 4 (bijrk+ bkjri- bikrj) + 3rirjrk]
I
C3 = - - - - (8.60)
4n(l-v)
C 4 = (1- 2v).
The tractions on a surface through x with outward normal ni can be calculated
from:
Tij (x, y) = Dikj (x, y) nk (x) (8.61)
which gives
C3
Tij (x, y) = 2R2 (C 4 (nj ri - ni rj) + [3 ri rj + C4 bij] r, n,) . (8.62)
x"u,
constant over a small area S. The displacements at any point x due to displace-
ments and tractions at surface S+ are given by:
(8.63)
where
L1 V ij (x, Yc) = S V ij (x, y) dS (y) (8.64)
and s
L1 Tij (x, yJ = S Tij (x, y) dS (y) (8.65)
s
At the dislocation we have,
tf = - ti (8.66)
Also, since S+ and S_ are at identical location but have opposite outward normals
L1 Vu = L1 Vi} and L1 TU = - L1 Ti} (8.67)
Therefore the displacements at x due to the displacements and tractions at both S+
and S_ are
Ui(X) = L1 Vij tt + L1 Uij r; + L1 Tij ut + L1 Tij uj (8.68)
Substitution of (8.67) gives:
Ui(X) = -L1Tij (ut- un =L1Tij dj (8.69)
where
dj = ut- uj (8.70)
Using equation (8.58) the stresses are given by
O'idx) = L1 Sikj dk (8.71 )
where
Cs
Sijk = -,3 {3 [C4 [)ij'k + V ([)ik'j + [)jk';)
.
- 5'i'j 'k] n",
Diering [35] has shown that the solutions obtained here are in fact identical to the
ones by Starfield [5). The derivations presented here assume that the displacement
discontinuity is parallel to the Xl - X2 plane. If this is not the case, the displace-
ments and tractions have to be transformed into local directions. Details are given
by Diering [35).
APPLICATIONS IN MINING 201
References
I Salamon, M D. G., Elastic analysis of displacements and stresses induced by the Mining
of seam or reef deposits. Part I Jnl. of S. Afr. Inst. Min and Metall, 64 (4), pp. 129-149,
1963
2 Salamon, M. D. G., Elastic analysis of displacements and stresses induced by the mining
of seam or reef deposits. Part II, Jnl. S. Afr. Inst. Min and Metall, 64 (6), pp. 197-218,
1964
3 Salamon, M D. G., Elastic analysis of displacements and stresses induced by the mining
of seam or reef deposits. Part III, Jnl. S. Afr. Inst. Min and Metall, 64 (10), pp. 468 - 500,
1964
4 Salamon, M. D. G., Analysis of displacements and stresses induced by the mining of seam
of reef deposits. Part IV, Jnl. S. Afr. Inst. Min and Metall, 65 (5), pp. 319- 338, 1964
5 Starfield, A M., Crouch, S. L., Elastic analysis of single seam extraction. In: New Horizons
in Rock Mechanics. (Eds. Hardy, H R., Jr. and Stefhanko, R.) ASCE, N.Y., pp. 421-439,
1973
6 Crouch, S. L., Solution of plane elasticity problems by the displacement discontinuity
method. Int. J. Numer. Meth. in Eng., 10, pp. 301- 343, 1976
7 Crouch, S. L., Computer simulation of mining in faulted ground. Jnl. S. Afr. Inst. Min and
Metall, 79(6),pp.159-173, 1919
8 Diering, T. A C, Simulation of mining in non-homogeneous ground using the displace-
ment discontinuity method. Jnl. S. Afr. Inst. Min and Metall, 80, No.7, pp. 169-173,
1980
9 Wiles, T. D., Curran, 1. H, A general 3-D displacement discontinuity method. Proc. 4th
Int. Con/. on Numerical Methods in Geomechanics, ed. Eisenstein, A A Blakema,
Rotterdam, 1, pp. 103-109, 1982
10 Crawford, A M, Curran, T. H, Higher order functional variation displacement dis-
continuity method. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 19, pp. 143-148,
1982
II Deist, F. H et aI., A new digital method for three dimensional stress analysis in elastic
media. Rock Mechanics 5, 189 - 202, 1973
12 Deist, F. H, Georgiadis, E., A computer system for three-dimensional elastic analysis
using a boundary element approach. Chamber of Mines of South Africa, Research Report
No. 43176, Project No. GS 1510
13 Rizzo, F. T., An integral equation approach to boundary value problems of classical
elastostatics. Q. Appl. Math. 25,83-95,1967
14 Cruse, T. A, Application of the boundary integral equation method to three-dimensional
stress analysis. Computer & Structures, 3, pp. 509- 527, 1973
15 Lachat, 1. C, A further development of the boundary integral technique for elastostatics.
Ph.D. Thesis, University of Southampton, u.K.
16 Lachat, J. C, Watson, 1. 0., Effective numerical treatment of boundary integral equations:
A formulation for three-dimensional elastostatics. Int. Jnl. for Num. Meth. in Eng. 10, pp.
991-1005,1976
17 Bannerjee, P. K., Butterfield, R., Boundary element method in Geomechanics. In: Finite
Elements in Geomechanics. Ed. Gudekus, G., Chichester, Wiley 1977
18 Brady, B. H G., Bray, 1. W., The boundary element method for elastic analysis of tabular
orebody extraction, assuming complete plane strain. Int. 1. Rock Mech. Min. Sci. &
Geomech. Abstr., 15,pp. 29-37,1978
19 Brady, B. H G., A boundary element method for three-dimensional elastic analysis of
tabular orebody extraction. Proc. 19th Rock Mech. Symp., pp. 431- 438, 1978
20 Cruse, T. A, Formulation of boundary integral equations for three-dimensional elasto-
plastic flow. Int. 1. Solids & Structures, 7
21 Hocking, G., Stress analysis of underground excavations incorporating slip and separation
along discontinuities. In: Recent Advances in Boundary Element Methods. (Ed. Brebbia,
C A) Pentech Press, Plymouth, pp. 195-214, 1978
202 APPLICATIONS IN MINING
22 Wardle, L. 1., Crotty, 1. M., Two Dimensional boundary integral equation analysis for
non-homogeneous mining applications. In: Recent Advances in Boundary Element Methods.
Ed. Brebbia, C. A., Pentech Press, Plymouth, pp. 233 - 251, 1978
23 Telles, 1. C. F., Brebbia, C. A., Plasticity. In: Progress in Boundary Element Methods. Ed.
Brebbia, C. A., Pen tech Press, London 1981
24 Beer, G. et aI., Efficient analysis in geomechanics. 4th lnt. Conf. on Numerical Methods in
Geomech. Ed. Eisenstein, A. A. Blakema, Rotterdam, 1, pp. 5-13, 1982
25 Beer, G., Meek, 1. L., Coupled Finite Element - Boundary Element analysis of infinite
domain problems. lnt. Conf. on Num. Meth. for Coupled Problems, ed. Bettess et aI.,
Pineridge Press, Swansea, u.K., 1981
26 Sinha, K. P., Displacement discontinuity technique for analysis stresses and displacements
due to mining in seam deposits. Ph.D. dissertation, University of Minnesota, 1979
27 Diering, T. A. c., An improved method for the determination by a MINSIM type of
analysis of stresses and displacements around tabular excavations. lnl. S. Afr. lnst. Min
Metall, 80, No. 12, pp. 225 - 228, 1980
28 Wilson, R. B., Cruse, T. A., Efficient implementation of anisotropic three-dimensional
boundary integral equation stress analysis. lnt. 1. Num. Meth. in Engng., 12, 1383 - 1397
29 Banerjee, P. K., Integral equation methods for analysis of piecewise non-homogeneous
three-dimensional elastic solids of arbitrary shape. lnt. 1. Mechanical Science, 18, pp.
293-303,1976
30 Brady, B. H. G., Bray, 1. W., The boundary element method for determining stresses and
displacements around long openings in a triaxial stress field. lnt. 1. Rock Mech. & Min.
Sci., 15 (I), pp. 21-28,1978
31 Hoek, E., Brown, E. T., Underground excavations in rock. lnt. of Min. and Met., 44 Port-
land Place, London WIN 4BR, England
32 Massonet, C. E., Numerical use of integral procedures. Chapter 10 in Stress Analysis. Ed.
Zienkiewicz, John Wiley, London 1965
33 Oliviera, E. R. A., Plane stress analysis by a general integral method. 1. ASCE, 94, No.
EMI, PP. 79-101,1968
34 Watson, J. 0., Advanced implementation of the boundary element method for two and
three dimensional elastostatics. Chapter 3 in Developments in Boundary Element meth-
ods-i. Ed. Banerjee, Applied Science Publ., 1979
35 Diering, T. A. C., Further developments of the boundary element method with applica-
tions in mining. MSc. thesis, University of the Witwatersrand, Johannisburg, South
Africa, 1981
36 Riccardella, P. c., An improved implementation of the boundary integral techniques for
two-dimensional elasticity problems. Carnegie lnst. of Tech., Dept. Mech. Eng., Pitts-
burgh, Report No. SM-72-76, 1972
37 Riccardella, P. c., An implementation of the boundary integral technique for planar
problems of elasticity and elasto-plasticity. Carnegie lnst. of Tech., Dept. Mech. Engng.,
Pittsburgh, Report No. SM-73-10, 1973
38 Crotty, J. M., Wardle, L. 1., Program BITEM user's manual. Institute of earth resources,
CSIRO Division of Applied Geomechanics, P.O. Box 54, Mount Waverley, Victoria,
Australia, 1980
39 Wardle, L. 1., Three dimensional boundary element program for mining applications:
MIN3D2, Technical Report No. 116, CSIRO Div. of Applied Geomech., Australia, 1980
40 Watson, J. 0., The solution of boundary integral equations of three-dimensional
elastostatics for infinite regions. Proc. First lnt. Seminar on Recent Advances in Boundary
Element Methods, Univ. of Southampton, 1979
41 Lachat, 1. c., Watson, 1. 0., Progress in the use of boundary integral equations, illustrated
by examples. Computer Methods in Applied Mechanics and Engineering, 10 (3), pp.
273- 289, 1977
42 Watson, J. 0., Hermitian cubic boundary elements for plane problems of fracture
mechanics, Res. Mechanica, 1981
43 Zienkiewicz, O. c., The Finite Element Method - Third Edition, McGraw-Hill 1977
44 Crotty, J. M., A block equation solver for large unsymmetric matrices arising in the
Boundary Integral Equation Method. Int. 1. Numer. Meth. in Engng., 18, pp. 997-1017,
1982
APPLICATIONS IN MINING 203
45 Hocking, G., Stress analysis of the 1100 orebody at Mount Isa. Proc. 19th U.S. Symposium
on Rock Mechanics, Stateline, Nevada, 1978
46 Crotty, 1. M., Wardle, L. T., Integral equation analysis of piecewise homogeneous media
with structural discontinuities. In preparation. C.S.I.R.O., Div. of Applied Geomechanics,
Australia, 1982
47 Bolteus, L., Tullberg, 0., BEMSTAT - A new type of boundary element program for two-
dimensional elasticity problems. 3rd Int. Sem. on Recent Advances in Boundary Element
Methods, Irvine, USA, 1981
48 Zienkiewicz, O. c., Pande, G. N., Time .dependent multi-laminate model of rocks - A
numerical study of deformation and failure of rock masses. Int. Jnl. N urn. Anal. Meth.
Geomech., 1, pp. 219- 247, 1977
49 Venturini, U. S., Brebbia, C. A., Boundary Element formulation for non-linear applica-
tions in geomechanics. 4th Int. Conf. on Recent Advances in Boundary Element Methods,
ed. Brebbia, Southampton, England, 1982
50 Venturini, u., Brebbia, C. A., The Boundary Element Method for the solution of no-tension
materials in Boundary Element Methods. Ed. Brebbia, C. A., Springer Verlag, Berlin, New
York 1981
51 Zienkiewicz, O. c., Kelly, D. W., Bettess, P., The coupling of the finite element method
and boundary solution methods. Int. J. Numer. Meth. Eng., 11 (2), pp. 355- 375
52 Mustoe, G. G. W., Symmetric variational boundary integral and finite element pro-
cedures in continuum mechanics. Ph.D. dissertation, Dept. of Civil Eng., Univ. of
Wales, Swanesa, U.K, 1980
53 Beer, G., Finite element, boundary element and coupled analysis of unbounded problems
in elastostatics. Int. 1. Numer. Meth. Engng., in press
54 Brady, B. H. G., Wassyng, A., A coupled finite element - boundary element method of
stress analysis. Int. 1. Rock Mech. Min. Sci. & Geomech. Abstr., 18, pp. 475 - 485, 1981
55 Beer, G., BEFE - Users manual. Dept. of Civil Eng., University of Queensland, St. Lucia,
Australia, 1982
56 Beer, G. et aI., Prediction of the behaviour of shale hanging walls in deep underground
excavations. 5th ISRM congress, Melbourne, Australia, 1983
57 Bettes, J. A., Economical solution technique for Boundary Integral matrices. Report
C/R/399/81, Dept. of Civil Eng., University College of Swansea, U.K
58 Brebbia, C. A., Walker, S., Boundary Element techniques in engineering. Newness-Butter-
worths, 1980
59 Mindlin, P. D., Force at a point in the interior of a semi-infinite solid. Physics, 7, pp.
195-202,1936
60 Telles, 1. C. F., Brebbia, C. A., Boundary element solution for half-plane problems. Int. J.
Solids Structures, 17, No. 12, pp. 1149-1158, 1981
61 Wardle, L. 1., Three-dimensional solutions for displacement discontinuities in cross-aniso-
tropic media. CSIRO, Div. of Appl. Geomech., Technical paper No. 34, 1980
62 Wardle, L. J., Fundamental solutions and boundary integral formulation for tabular
excavations in a non-homogeneous rock mass. CSIRO, Div. of Appl. Geomech., Techn.
report No. 133,1982
63 Wardle, L. 1., Boundary element methods for solution of stress analysis problems in
mining. CSIRO, Div. of Appl. Geomech., Techn. Report No. 102, 1980
64 Meek, 1. L., Matrix Structural Analysis. McGraw-Hill Co, New York 1971
65 Hebblewhite, B. K, The use of modelling techniques for mine design. Australian Journal
of Coal Mining Technology and Research, No. I, 1982
66 Brebbia, C. A., Georgiou, P., Combination of boundary and finite elements in elasto-
statics. Appl. Math. Modelling, 3, pp. 212 - 220, 1979
67 Banerjee, P. K, Butterfield, Boundary Element Methods in Engineering'Science. McGraw-
Hill, London 1981
68 Bywater, S., Cowling, R., Black, B. N., Stress Measurement and Analysis for Mine Plan-
ning. 5th ISRM Congress, Melbourne, Australia, 1983
69 Beer, G., BEFE - a combined oundary element finite element computer program. Adv.
Eng. Software, Vol. 6, No.2
Chapter 9
9.1 Introduction
The Boundary Element Method (BEM) is now an effective tool for the numerical
analysis of nonlinear as well as linear problems. Some of those problems are
physical or material nonlinear, such as elastoplasticity and creep in solid mechan-
ics [1]. Another type of nonlinear problem are those concerned with geometrical
nonlinearities, such as finite deformations. Among a variety of geometrically
nonlinear behaviours in solid mechanics, the finite deflection of flat plates or shells
is one of the most important problems from the standpoint of engineering practice.
Very little work has been carried out on the treatment of geometrically nonlinear
equations using BEM. The conventional von Karman nonlinear differential
equations for finite deflections of plates represent a complex situation involving
coupling between inplane and out-of-plane deformations. So far the BEM appears
to have been applied to find approximate solutions for the finite deflection of
plates and shells [2-5). In these studies, a pseudo-linear governing equation known
as the Berger equation [6] has been used, which approximates nonlinear bending
behaviour without introducing the above-mentioned coupling. The Berger
decoupled field equation is not based on a sound physical assumption and is
occasionally referred to as the Berger "hypothesis". While the equation may not
theoretically be precise, several have shown its effectiveness for practical applica-
tions provided that the inplane displacements are constrained on the boundary.
The present chapter discusses the previous works and presents the general
equations for finite deflection of thin elastic, flat plates based on the Kirchhoff-
Love theory. The chapter also shows an integral equation formulation of the von
Karman-type nonlinear governing equations and a short derivation of the approxi-
mate governing equation due to Berger and its corresponding integral form. Some
results illustrate the application of the approximate integral formulation for the
Berger equation. Finally, the chapter deals with the extension of the BEM to non-
linear bending problem with thin shallow shells and sandwich plates and shells.
Some representative numerical results are included.
Boundary r n
Xl
Fig. 9.1. Flat plate and its element. Coordinate and notations
the Kirchhoff-Love hypothesis and referred to the Cartesian coordinates x), X2, Z
shown in Fig. 9.1. The summation convention is used throughout ranging from I
to 2. Linear and nonlinear terms are distinguished by the superfixes I and n respec-
tively (for the complete notations, see Appendix III).
I) Fundamental definitions for plates
(1) Membrane strains:
eij = efj + elj (9.1)
elj = t (UiJ + Uj,i) (9.2)
2 ,I ,J.
e~·=.!.w,w
IJ (9.3)
V = Qi ni = Vi + vn (9.16)
Vi = Mij,jni (9.17)
vn = N ij W,j ni (9.18)
Pi =Pi (i=1,2)
Mn = Mn (9.25)
Kn =Kn
(2) Geometrical boundary conditions on r 2:
Ui = Ui (i= 1,2)
W = W (9.26)
W,n= w,n
where overscored values are prescribed on the boundary.
FINIlE DEFLECTION OF PLATES 207
and
J[- DV 4 W + Q7 (Uk. w),;+ .0] w* dQ (9.35)
Q
The field variables other than the displacements and moments Mn can now be
defined and composed of linear and nonlinear parts,
N;j (Uk. w) = N;~ (Uk) + Nij (w) (9.36)
p; (Uk, w) = pi (Uk) + p7 (w) (9.37)
Kn (Uk. w) = K~ (w) + K~ (Uk. w) (9.38)
Q; (Uk. w) = QI (w) + Q7 (Uk. w) (9.39)
We can take the same weight functions uT for Eqs. (9.32) and (9.34), and w* for
Eqs. (9.33) and (9.35), respectively. Consequently, the related kernels are identical
to those used in the linear decoupled case,
pT (un = 1;* (un, Nfj (un = Nfl (uk), K: (w*) = K~* (w*) . (9.40)
We can now integrate Eqs. (9.34) and (9.35) by parts taking into consideration
Eqs. (9.36) to (9.40). This gives
and
- D JwV 4 w* dQ + J.ow* dQ - J[Mn(w) w~n- w,nM: (w*) + Kn (Uk. w) w*
Q Q r
- w K: (w*)] dr = JQ7 (Uk. w) w:'i dQ (9.42)
Q
Identical results could be obtained using weighted residuals for the governing
equations and boundary conditions as shown in Reference [13]. The main
difference between Eqs. (9.41)-(9.42) and Eqs. (9.32)-(9.33) consists of the
additional nonlinear terms written on the former right hand sides. These terms are
given by the integrals on the domain Q, with terms of win Eq. (9.41) and of wand
Uk in Eq. (9.42) representing the coupling between inplane and out-of-plane
deformations.
The weight functions introduced in Eqs. (9.32), (9.33) and (9.41), (9.42) are the
respective fundamental or principal solutions of the following equations:
Nt [u~ (P, Q)l,j + bm (P, Q) =0 (9.43)
and
V 4 w* (P, Q) + b (P, Q) = 0 (9.44)
These solutions are known as the Kelvin solution for the two-dimensional infinite
body in a state of plane stress and as the fundamental solution to the two-
dimensional biharmonic equations for the out-of-plane deformation, respectively.
Equation (9.43) represents the two-dimensional field at a point P under unit
inplane force concentrated at point Q acting in the m-direction (m = 1,2).
Equation (9.44) corresponds to the deflection at P due to a unit concentrated
lateral load at Q. Their solutions are well-known and functions of the distance r
FINITE DEFLECTION OF PLATES 209
I+v [ rn (I-V
p%[u%(P,Q)]m=-4- --(jkm+ 2r k r m) ---(rmnk-rknm)
I-v ] (9.47)
nr ' I+v " l+v' ,
1
w* (P, Q) = - ~ r210g r (9.48)
where Ukm and Pkm denote the components of u% and p% due to the unit inplane
force applied at Q in the m-direction. The related kernels of Eq. (9.42), w;n' M:
and K:,can be obtained by differentiating Eq. (9.48) (for further detail, see Ref.
[II]).
Substituting the corresponding fundamental solutions into Eqs. (9.41) and (9.42)
and computing the Cauchy principal value integrals, we finally have
- C (P) Ui(P) + S{pdum(Q)] u% (P, Q) i - udQ) p% [u: (P, Q)] J dT
r
= SN[m [w (Q)] u% (P, Q)i,m dQ (Q E r, Q E Q) (9.49)
Q
and
D C(P) w(P) + JP(Q) w* (P, Q) dQ
Q
where the value of the coefficient C (P) depends on whether the point P belongs to
the domain Q or to the smooth boundary r, i.e.,
C(P) = U (P
(P
E
E
Q)
T)
(9.51 )
where Ai, ... , Di , E i , ... , Hi (i = 1,2) are influence coefficient matrices and the
right hand side vectors Ni (i = 1, ... ,4) contain the nonlinear and inhomogeneous
terms. To compute the nonlinear terms an iterative procedure is required. The
(t + l)-st approximation is estimated using the t-th approximation for the right
hand side terms. Calculation can then be repeated until one reaches the required
accuracy. Equations (9.55) to (9.58) are required only on the boundary, i.e.,
C (P) = ~ in Eqs. (9.49) and (9.50). However, to calculate the right hand side terms
Ni one needs to compute the internal values of the displacements and their
derivatives for which Eqs. (9.49) and (9.50) are used with C (P) = I. These integra-
tions are usually performed numerically.
and quasi-linear it has been extensively applied in place of the complete equation
derived in section 3. The equation has been used for dynamic as well as static
problems, shallow shells, sandwich plates and several others [14-18]. It is impor-
tant to notice that the Berger solution can give a fairly good approximation to the
problem provided that the inplane displacements are restrained on the boundary.
If a thin elastic plate is SUbjected to a lateral load P (XI, X2) and a thermal field
T (XI, X2, z), the total potential energy can be expressed as
The complete von Karman-type equations can be derived from Eq. (9.59) using
variational calculus. If we neglect h, i.e., the second invariant of the membrane
strains, as suggested by Berger, and then apply variational calculus, the resulting
equation is,
(9.63)
11 -
(1 + v) NT = - -h=
')(22
const. (9.64)
Eh 12
Equation (9.63) is linear with respect to w, but the parameter ')(2, sometimes called
Berger constant, is dependent on the external force or temperature field. In other
words, the second term in Eq. (9.63) represents nonlinear relationship between
deformation and intensity of external actions.
The magnitude of the Berger constant ')(2 can be estimated as follows. Let us
integrate both sides of Eq. (9.64) over the whole two-dimensional plate domain Q,
i.e.,
x2=
12
h2Aok
[ 1
UI,I+U2,2+'2(W,I+W,2)-
2 2 (1+v)
Eh
NT] dQ. (9.66)
212 FINITE DEFLECTION OF PLATES
Since now the left hand side contains only the biharmonic operator in terms of the
deflection and a harmonic term in MT exists on the right hand side, the following
generalized Green identity will hold for smooth boundary [19]:
MT
MI(w) = Mn(w) - - - (9.70)
I- v
By using the known fundamental solution (9.48) and applying Cauchy principal
values, we arrive at the following equation:
D C(P) w(P) - J{w* (P, Q) K~T[w (Q)] - w;n (P, Q) MI[w (Q)]
r
+ W,n (Q) M: [w* (P, Q)] - w (Q) K~· [w* (P, Q)]} dr
(Q E r, Q E Q) (9.71)
where the coefficient C (P) depends as previously on the position of the point P.
Another equation for the normal derivative of the deflection on the boundary is
obtained by differentiating Eq. (9.71) in the no direction, i.e.,
tD w,no(P) - J{w~no(P, Q) K~T[W(Q)] - w~nno(P, Q) Mr[w(w(Q)]
r
- w,n(Q) ~ [w*(P, Q)),no- w(Q) K~· [w* (P, Q)),no} dr
t
= [{D x2 [V 2w(Q)] + p(Q)} w~no(P, Q) - l~Tv (Q) V2w~no(p, Q)] dQ
(P, Q E r, QE Q) (9.72)
FINITE DEFLECTION OF PLATES 213
x2 = ~ J
h 2AD D
[2..2 (2 + W,22) _ (1 +E v)h NT] dQ
W,I
(9.75)
x2--~J[2..WV2W+
2
(1 +V)NT]dQ (9.77)
- h AD D 2 Eh
In this section some numerical results obtained using the integral formulation of
Berger equation are presented.
214 FINITE DEFLECTION OF PLATES
The plate boundary and domain are discretized into constant elements and cells.
In this way a set of algebraic simultaneous equations is obtained from Eqs. (9.73)
and (9.74) using numerical integration. Calculation of the integral related to the
Berger constant can also be carried out numerically for a given distribution of the
deflection and thermal field. An iterative procedure is applied until the required
convergence is achieved.
To check the validity of the proposed approximate procedure, some simple
circular and square plate problems have been solved. For these problems, the
nonlinear solution can be obtained using conventional methods such as Runge-
Kutta-Gill (RKG) method or classical energy method due to their simple
geometry. The plates are considered to be subjected to uniformly distributed
lateral load
p = const
or a distributed thermal field given by (say, for circular plates),
l
T(R, z) = To + T] (I - ~:)]( 1+ ~ ~), R2 = XI + x~
where To and T] are constants. The temperature is assumed to vary linearly within
the thickness of the plate.
Figures 9.2 to 9.5 show results for the maximum deflection of circular and
square plates; obtained by using 24 boundary elements and 72 cells in the case of
circular plate, and 48 boundary elements and 64 cells for square plates [2, 3]. The
boundary conditions considered are indicated in each of the figure captions.
Poisson ratio has been taken v = 0.3.
The present boundary element solutions (marked by small circles in figures) are
compared against the more rigorous results obtained using the von Karman
equation. A set of results obtained using coarse discretization (12 boundary
elements) is also included in Fig. 9.3 and plotted as black small circles.
1.0 1.0,----------------,
Linear
Linear solution
0.8 0.8 conventional)
solution equation
(RKG Method)
t 0.6
~ '
~:ventional
equation
t 0.6 Berger
equation
(RKG method)
,JO.4 (energy method)
0.2 0.2
o Berger equation (BEM) • 0 Berger equation (BEM)
0 10 20 30 40 50 60 60
pa4fDh -
Fig. 9.2. Maximum deflections of immov- Fig. 9.3. Maximum deflections of immov-
ably clamped square plate under uniform ably clamped circular plate under uniform
lateral load lateral load
FINITE DEFLECTION OF PLATES 215
0.4,------------r-----.,
Linear
solution
0.3 Conventional
equation
(RKG method)
~
10.2 Berger
equation
(RKG method)
,J
0.1
4
pa 4l D h -
Fig. 9.4. Maximum deflections of immovably supported, rotation free circular plate under
uniform lateral load
Berger
equation
(RKG method)
o O.B 1.0
WmG,Ih -
Fig. 9.5. Maximum deflections of immovably clamped heated circular plate
The Berger solution for 24 boundary elements agrees well with the more
accurate results. A significant improvement in the accuracy of the solution is
noticeable when comparing this solution with the one obtained using the coarse
mesh (Fig. 9.3). It should be pointed out that it is always difficult to carry out
rigorous nonlinear analysis for plates with arbitrary plans, hence the importance of
boundary element formulations to solve this problem.
Figure 9.5 shows the deflections of a plate heated by inhomogeneous temper-
ature distribution. The boundary element solution coincides with the RKG
solution of the Berger equation for the whole range of the temperature parameter
TJITo. The difference between the results obtained using the Berger or the
complete von Karman formulations, i.e., the uncoupled and coupled solutions,
varies with the magnitude of the parameter TJITo and the stage of deformation.
Nevertheless, the difference does not seem significant from the practical point of
view and the Berger solution may trace well the coupled results even for highly
nonlinear behaviour.
216 FINITE DEFLECTION OF PLATES
There have been many applications of Berger original idea to approximate the
finite deflection of shallow shells and sandwich structures. In this section, we will
outline briefly how this can be extended to find an integral formulation for
problems with shallow shells and sandwich plates and shells.
Boundaryr S n
Xl
kl/kl
dXlydXl
~
Fig. 9.6. Shallow shell and its element. Coordinate and notations
FINI1E DEFLECTION OF PLA1ES 217
The final integral formulae are of the same form as those in Eqs. (9.71) and (9.72)
except that the term V 2 w has to be replaced by V 2 w + k] + k 2 • It must be pointed
out that the Berger constant for shallow shells is also affected by the curvature.
The results obtained for axisymmetric finite deflections of circular spherical
caps under constant lateral load and distributed temperature are shown in Figs. 9.7
and 9.8. Berger results coincide with the corresponding RKG solutions for a whole
range of total curvature.
1.2.------------.---------,
1.0
Linear
solution
0.8 (flat plate)
-0.2
0.4
-0.4
0.2
o 4 6 8 10
pa"lEh 4- -
Fig. 9.7. Maximum deflections of immovably clamped spherical cap under uniform lateral
load
2.5 r - - - - - - - - - - - - . - - - - - - - - - - - - - - ,
a 1 +k2) =0.4
2
7I(k
~1.0 0.2
o
-0.2
0.5 -0.4 o BEM solution
- RKG solution
OL-_~_~_~_~_~ _ _L_~_~_~~~
z
Fig. 9.9. Sandwich shell and its element. Coordinate and notations
Since the principal term of Eq. (9.81) is the biharmonic representation of the
deflection, we may follow a procedure analogous to that of section 7.1, to obtain
the following integral equations:
C(P) w(P) + f{K~ (Q) W* (P, Q) - M~(Q) w;n (P, Q) + ~ [w* (P, Q)] w, n(Q)
r .
-K~'[W*(P,Q)]W(A)}dr=-s[ A x~2+ I 12
(['1 2 w*(P,Q)]w(w(Q)
=- b[A x~+ I (['1 2w;no(P, Q)] W(Q) + w;no (P, Q)[k (Q) + k2 ({2)]} I
where
C(P) = n (P E Q)
(P E r)
(9.88)
(9.89)
From these formulae the discretized algebraic equations can be obtained as usual
and a numerical iterative scheme can be carried out in order to obtain the solution.
Some numerical results are shown in Figs. 9.10 to 9.13 for clamped circular flat
sandwich plates and square shallow sandwich shells under uniform lateral load and
specified temperature difference between both faces. The following geometrical
and material parameters are used in the calculation:
t = 0.635mm, h = 17.13 mm, vI = 0.3
EI = 0.72 x 10 5 MPa, G8 = 0.414 x 102 MPa
Temperature distribution for circular sandwich plates:
.10-3
1.5.---------------r-. 6
5
1.2
J O.6
0.3
o BEM solution o BEM solution
- - RKG solution
1.5 .--------~---~
p
~
.::::;,::::::::::,.::::{
1.2 1---20--1
/
Linear
solution
(flat plate)
4a 2
T(k,+k2 )=
rrh
0.04
0.02
o
OJ -0.02
-0.04
~2a~
Linear
solution
(flat plate)
~ 1.':-0---::'::----:-':----:"----:L-----'----:L----,l-----'----:L-----'1.0
WmQ/h -
Fig. 9.13. Maximum deflections of heated shallow square sandwich shell (Gc = G8)
In this chapter, the application of BEM to solve finite deflection of elastic thin
plates, shallow shells and sandwich plates and shells has been presented. This is an
important problem in structural and solid mechanics for which integral equation
formulation can be conveniently applied.
Integral equation formulations based on the von Karman theory have been
described as well as the so-called Berger approximation. The former gives a
coupled system of field equations which is reduced to an uncoupled system for the
latter.
Some numerical solutions are presented for the equation under the Berger
hypothesis. This approximation gives fairly good results provided that the
boundary conditions are suitable and represents a simple way of simulating the
complex nonlinear behaviours of plates and shells.
Although the results presented in this chapter were restricted to simple
geometrical shapes of plates and shells, external force distribution, special
boundary condition - thus for the case of Berger equation - piecewise constant
elements and smooth boundaries; the potentialities of the method are demon-
strated. Further studies should be carried out to implement the technique in a
more general form. Moreover, other integral formulations could be studied in
addition to those given by the Berger and von Karman field equations.
Acknowledgement
The authors wish to express their cordial thanks to Dr. C. A. Brebbia of
Southampton University and Institute of Computational Mechanics for his
valuable discussions and suggestions.
222 FINITE DEFLECTION OF PLATES
References
I Brebbia, C. A, Progress in Boundary Element Methods. Pentech Press, London 1981
2 Kamiya, N., Sawaki, Y., An integral equation approach to finite deflection of elastic
plates. Int. J. Non-Linear Mech., 17,pp.187-194, 1982
3 Kamiya, N., Sawaki, Y., Nakamura, Y., Fukui, A, An approximate finite deflection
analysis of a heated elastic plate by the boundary element method. Appl. Math. Modell.,
6,pp. 23-27,1982
4 Kamiya, N., Sawaki, Y., A simplified nonlinear bending analysis of flat plates and shallow
shells by boundary element approach based on Berger equation. Proc. Int. Conf. Num.
Meth., pp. 289- 297, 1982
5 Kamiya, N., Sawaki, Y., Nakamura, Y., Boundary element nonlinear bending analysis
of clamped sandwich plates and shells. Proc. Fourth Int. Conf. on Boundary Element
Methods in Engineering, Brebbia, C. A, ed., pp. 515-525, 1982
6 Berger, H. M, A new approach to the analysis of large deflections of plates. J. Appl.
Mech., 22,pp. 465-472,1955 .
7 Timoshenko, S., Woinowsky-Krieger, S., Theory of Plates and Shells. McGraw-Hill, New
York 1959
8 Washizu, K, Variational Methods in Elasticity and Plasticity. Pergamon Press, 1968
9 Kamiya, N., Sawaki, Y., Integral formulation for nonlinear bending of plates - Formula-
tion by weighted residual method. Zeit. ang. Math. Mech., 62, pp. 651-655,1982
10 Bergman, S., Schiffer, M., Kernel Functions and Elliptic Differential Equations in Mathe-
matical Physics. Academic Press, New York 1953
II Bezine, G. P., Gamby, D. A, A new integral equation formulation for plate bending
problems. Recent Advances in Boundary Element Methods. Brebbia, C. A, ed., Pentech
Press, pp. 327 -342, London 1978
12 Stem, M, A general boundary integral formulation for the numerical solution of plate
bending problems. Int. J. Solids Struct., 15, pp. 769 -782, 1979
13 Brebbia, C. A, The Boundary Element Methodfor Engineers. Pentech Press, London 1978
14 Nash, W. A, Modeer, J. R., Certain approximate analyses of the nonlinear behaviours of
plates and shallow shells. Proc. Symp. Theory of Thin Elastic Shells (IUTAM), pp. 331-
354, 1960
15 Nowinski, J. L., Ismail, I. A, Certain approximate analyses of large deflections of
cylindrical shells. Zeit. ang. Math. Phys., 15, pp. 449-456, 1964
16 Pal, M C., Large deflections of heated circular plates. Acta Mech., 8, pp. 82-103, 1969
17 Kamiya, N., Governing equations for large deflections of sandwich plates. AIAA J., 14,
pp. 250-253,1976
18 Kamiya, N., Analysis of the large thermal bending of sandwich plates by a modified
Berger method. J. Strain Anal., 13, pp. 17 - 22, 1978.
19 Kamiya, N., Sawaki, Y., Nakamura, Y., Thermal bending analysis by boundary integral
equation method. Mech. Res. Comm., 8, pp. 369-373, 1981
Appendices
Subscripts
i, j, k ... (except n, s) I, 2 (summation convention for repeated indices)
(),i 0 ()/OXi
(),n o( )/on normal derivative on the boundary
(),s o( )/os tangential derivative on the boundary
Superscripts
b refer to bottom face
c refer to core
f refer to face
I refer to linear term
m refer to average value on top and bottom faces
n refer to nonlinear term
refer to top face
* weight function or its related functions
Chapter 10
Trefftz Method
by Ismael Herrera
10.1 Introduction
10.2 Scope
aQ
Fig. 10.1
IREFFIZ ME1HOD 227
where Yo stands for the trace of U on oQ. In general, for simplicity the symbol Yo
will be omitted when it is clear from the context that we refer to the boundary
values. It can be noticed that the linear space D defined by (IO.5) is not closed.
Indeed, a metric is not defined in the whole space.
Let
N p = {u E DiLl u = 0 in Q} (10.6)
and
1= {u E Diu = 0 on 0 Q} (10.7)
The first of equations (l0.8) is equivalent to (10.3 a), while the second one to
(l0.3b).
A first advantage of formulating the problem in this manner is connected with
its existence properties. Clearly equation (l0.3 b) is equivalent to u - U = V - U on
oQ. By well known results on the existence of solution [41], this problem possesses
a unique solution. Indeed, given U E D and V E D there are real numbers rand s
such that U E H r (Q) and the trace Yo (V - U) E H S (oQ). Then, u - U E H S + 112 (Q).
Therefore, u = U + (u - U) belongs to HI (Q) where t = min {r, s + 1I2}. This shows
u E D.
The above discussion also shows that there is no lack of generality by restricting
attention to the homogeneous case; i.e.
Llu = 0; on Q (l0.9a)
and
u = fag; on oQ (l0.9b)
The boundary method to be applied depends on the continuity of the solutions on
their boundary values. In principle it can be applied when the space of admissible
functions D is given by (l0.5). However, this would lead to consider inner products
in the space of boundary values H S (oQ) with arbitrary s, which may be
inconvenient in numerical applications. It is preferable to keep the computations in
../2(oQ) = HO(oQ), which, as will be seen, leads to least-squares fitting. This, can
be achieved if attention is restricted to functions with boundary values belonging
to HO(oQ) =../2 (oQ). When this condition is incorporated in the definition of the
space of admissible functions, one gets
(10.10)
This is again a linear space which is not closed.
In addition, in many applications it is necessary to compute the normal
derivative oulOn on the boundary oQ. Similar considerations lead to require that
ou/on belong to HO(oQ) =..f2(oQ). When these two requirements are incorporated
in the definition of the space of admissible functions, equation (10.5) becomes
(lO.Ila)
Here, as it is costumary, Yl u stands for the trace of the normal derivative on oQ.
This is again a linear space.
General results on the existence and continuity properties of solutions of elliptic
equations [41], imply that any harmonic function u whose trace Yo u belongs to
HO(oQ), necessarily is a member of HII2(Q). Therefore Npc HI/2(Q) in this case.
Even more, due to the continuity properties just mentioned, N p is a closed
subspace of HII2 (Q). This will be represented by N I12 (Q). Thus
Np=NI/2(Q) (lO.llb)
when D is defined by equation (10.10). Similarly, when equation (10.11 a) holds
corresponding properties imply that
Np=N3/2(Q) (lO.llc)
where N 3/2 (Q) is the subspace of harmonic functions belonging to H 3/2 (Q) which
can also be shown to be closed.
1REFFIZ METHOD 229
If qJ= {WI. W2, ... } C N 312 is a system of harmonic functions in Q which spans
N /2 then it also spans N 1I2 , because N 312 c N" 2 is dense in N" 2. Thus, let
3
qJ= {WI, W2, ... } C N 312 be such a system. Then, given U E N" 2 there is a sequence
of approximations
N
UN = L. a~ W n ; N = 1,2, ... (10.12)
n=1
such that
(10.l3)
Notice that ~ depends on the number of terms N, of the approximation. This is
essential in order for (10.13) to hold. When at; is independent of N, (10.12)
becomes a series and the approximation by a series can only be granted when the
system of functions qJ is orthogonal. This fact explains some of the difficulties that
were encountered in applications to electromagnetic field studies [43].
In order for representation (10.12) to be useful, it will be required to have a
procedure for deriving the coefficients at; from boundary data only. This is,
indeed, possible. General results on the existence and continuity properties of
solutions of elliptic equations [41] show that the range of the traces Yo u when
u E NII2(oQ) is Ho(oQ). Given u E N"2(Q) if the coefficients at; are chosen so
that
(10.14)
then (10.13) holds necessarily.
If qJ= {WI, W2, ... } c N3/2(Q) spans N312(Q), then the continuity properties
of elliptic equations imply that the associated system of traces {Yo WI, Yo W2, ... }
span HI (oQ), which in tum implies that
{YOWI,YOW2, ... } spans HO(oQ) (10.15)
These remarks show that the coefficients a~ can be chosen so that (10.14) holds.
Clearly, to this end it will be sufficient to take Yo UN as the projection of the
boundary values Yo u E HO (oQ) on the subspace spanned by {WI, ... , WN}. There-
fore, the coefficients can be computed by standard procedures for projecting on a
subspace. This yields a least-squares method because only ../2 (oQ) = HO (oQ) inner
products are being used.
This procedure leads to the system of equations:
N
In order for this system to be invertible it is required that the system of traces
{WI. ... , WN} C ~(oQ) be linearly independent. Then, the matrix Mnm can be seen
to be Hermitian.
Computation of the boundary values may need a special device. If the normal
derivatives ou/on are required on the boundary, in the case of Dirichlet problem,
230 TREFFIZ METIIOD
they can not be obtained from the approximating sequence UN, directly. Indeed, as
mentioned, from UN -> fao on HO (oQ), one can only grant that UN -> u in HI/2 (Q).
ou N ou ou N
Hence - - -> - in the metric of H- I (oQ) [41]; thus, - - diverges in
on on ou on
HO (oQ) in general. As a matter of fact, a;;
E H- I (oQ) is only defined as a con-
where we have written (,)a for the inner product in HO (oQ). Also <A u, v> = 0
whenever u E N p and v EN p. This yields the reciprocity relation
When u E N 3/2 (Q) the trace YI u (i.e. the normal derivative on the boundary)
spans
(10.22)
Here, {l}-L stands for the orthogonal complement in HO (oQ) of the constant func-
. ou . . ou
bon. Thus, -0 E {l }-L, If and only If, - E HO (oQ) and
n on
ou
S -0 dx =
ao n
o. (10.23)
Taking !!1J= {WI. W2, ... } c N 3/2 (Q) as before, and in view of (10.15), given any
u E N312(Q), it is possible to construct an approximation satisfying (10.21). Again,
for such representation to be useful it will be required to have a procedure for
computing b~ using boundary data only. This will be based on the use of the
reciprocity relation (10.20). It has interest to observe a general feature of re-
presentation (10.21); namely, that the normal derivative ouliJn of the sought
solution is not approximated by the normal derivatives of the basic system f%J, but
TREFFIZ METHOD 231
(10.24)
is minimized.
This leads to take the projection of ou/ on on the space spanned by
{WI, ... , WN} c Ho(oQ). This requires the orthogonality condition
ou .; N
( ;;--L.bnwn,Wm)
=0, m= 1, ... ,N (10.25)
un n=1
to be satisfied. Expanding (10.25), one gets
N
L Knm b~ = dm (10.26)
n=1
where
Knm= S wnwmdx; n,m= 1, ... ,N (10.27 a)
oQ
and
OU oW m
dm= S -0 wmdx= SfoQ-o-dx; m= 1, ... ,N (10.27 b)
oQn oQ n
Observe that the use of the reciprocity relation (10.20) has permitted to express dm
in terms of boundary data only.
An additional point must be mentioned. In order for the approximating
N OU
sequence L b~ Wn to be convergent, it is necessary that the solution -0 E HO (oQ).
n=1 n
This is granted if foQ E HI (oQ). Alternatively, this condition can be expressed in
matrix form. Let KN be the NxN square matrix whose elements are given by
(10.27 a). Similarly dN is the lxN vector defined by (10.27 b). Assume, the system
of traces {WI, ... , WN} C HI (oQ) is linearly independent, which is required in order
for the system (10.26) to be invertible, and denote by (KN)-l the inverse of KN.
Then, the sequence of real numbers
(10.28)
L d~< 00 (10.29)
n=l
The treatment of Neuman problem is similar. Let the space of admissible func-
tions be given again by equations (10.11 a). Then equation (I 0.9 b) is replaced by
OU
a;; = goQ; on oQ (10.30)
232 TREFFIZ METHOD
where the boundary values goQ E {l}.lc HO(oQ). The previous argument still
holds if (10.17) is replaced by
OWn OW m
S -;)-
Mnm = -;)- dx (10.31 a)
oQ un un
and
oWm
Cm = S gOQ-;)-dx (10.31 b)
oQ un
(10.34)
OWn OW m
S -;)-
Knm = -;)- dx (10.35 a)
oQ un un
and
dm =
oWm
S u -;)-dx = S goQ Wm dx (10.35 b)
oQ un oQ
Again, equations (10.26) have to be satisfied. When this is the case the solution u
in (10.34) fulfills (10.33). This method can be used to accelerate the convergence of
the approximating sequence on the boundary. As a matter of fact, when the system
Generally, when dealing with partial differential equations only some boundary
values of the functions and their derivatives are relevant in the discussion of the
problems. For example, for Laplace equation these are the function u and its
normal derivative au/an. For Elasticity the displacements u and tractions T (u).
When a boundary value problem is formulated, only one part of this boundary
information is prescribed and the other part must be derived after the solution has
been obtained. For Dirichlet problem, for example, u is prescribed and au/an is
derived. The converse has to be done in the case of Neuman problem. Approxi-
mating sequences for the complementary boundary values which depend on reci-
procity relations, such as (10.20) can be derived for very general classes of
differential equations. The reciprocity relations can be obtained from corre-
sponding Green's formulas. For example, from
au av }
f{ v Au - u A v} dx = S {v -a - u -a dx (10.36)
D aD n n
one obtains au av
S v -a dx = S u -a dx (10.37)
aD n aD n
Denote by .ff c *' the image of N 312 (Q) under the mapping u --+ it =
[Yo u, y, u] E*- It can be shown that .ff c *'
is closed in the metric of *- Notice
that the reciprocity relation (10.20) becomes (we assume Hilbert-spaces are being
taken with real coefficients):
(VI. U2) = (V2' u,) Y it E.ff &VE.ff (10.50)
A system of function f!Ij c N 3/ 2 (Q) will be said to be T-complete * if, for every
it E *. one has
(10.51)
Using this notation the following characterization of complete systems holds [19,
34].
Theorem 10.1. Let f!Ij c N 312 (Q). Then the following assertions are equivalent:
(i) f!ljc N312(Q) spans N312(Q) in the metric of H312(Q);
(ii) ;j c*' spans.ff in the metric of*'·
(iii) ;j c.ff is a T-complete system;
(iv) Equations (10.41) are satisfied when the spans are taken in the../ 2 (oQ) sense.
An advantage of having a system which satisfies any of the criteria (i) to (iv), is
that the same system can be used for both a Dirichlet and a Neuman problem.
Indeed, the same T-complete system can be used for any linear boundary condition
which is prescribed point-wise. Such condition can be written as
(10.52)
The arguments presented previously, can be extended to this case by introduction
of more general Green's formulas. This will be discussed in Section 3.
It has interest to observe that it is possible to develop systems which are
complete in regions which are, to a large extent, arbitrary. For example, the system
of harmonic polynomials given by (10.38) and (10.39), IS T-complete in any
bounded and simply connected region [29]. Also, the system
{Logr, Rez- n, Imz- n ; n = 1,2, ... } (10.53)
is T-complete in the exterior of any simply connected and bounded region which
contains the origin.
To develop general criteria establishing conditions under which a system which
is complete in a region is also complete in another one, is quite valuable.
Especially if such criteria are applicable to a wide class of partial differential
equations. For this purpose the notion of T-completeness is useful.
'" Trefftz-complete. Previously, such systems had been called c-complete by the author.
236 TREFFIZ METHOD
theory which permits obtaining such formulas systematically and which in some
respects enlarges the kind of problems that can be treated in this manner, has been
developed by the author [19, 28, 34]. The fundamental notions are closely related
with simplectic geometry [44].
Basically, what is done is to characterize the space of boundary values which are
relevant for each differential equation or system of such equations. Then such
space is decomposed into two subspaces. With every Green's formula there is
associated such a decomposition and conversely with every decomposition there is
a unique Green's formula. A procedure for reconstructing the Green's formula
when the decomposition is known, is established [34].
We consider a bilinear functional P defined on an arbitrary linear space D; it
will be denoted by P: D ~ D* because it can be thought as an operator defined on
the linear space D and taking values on its algebraic dual D* (this is the space of
linear functionals defined on D) [45]. The value of such bilinear functional at
elements U E D and v E D, will be denoted by (P u, v). The transposed bilinear
functional of P: D ~ D*, will be P* : D ~ D*; thus
(p* u, v) = (P v, u) (10.54)
The theory is applicable to general non-symmetric linear operators, although its
application to formally symmetric ones is simpler, because it does not require the
introduction of a formal adjoint. Here, attention is restricted to such operators.
Given an operator P: D ~ D*, we define the antisymmetric bilinear form
A= P- p* (10.55)
The operator A, given by (10.55) plays a central role in the theory. Firstly, we are
going to use it, to define the relevant boundary values. For this purpose, we
consider the null subspace N A of A; i.e.
(10.56)
With reference to the reduced wave equation
Au + k2U = 0 , on R (10.57)
as an example (recall that Laplace equation corresponds to the case k = 0),
consider the bilinear functional P: D ~ D*, given by
(P u, v) = S v (Au + k 2 u) dx (10.58)
R
Then A = P - p* is
(A u, v) =
ou - u a;;
lR { a;;
V
ov} dx (10.59)
The null subspace N p, is the linear subspace of functions which satisfy (10.57).
There are many ways of taking the linear space D. A convenient one is by means
of equation (10.11 a). This defines a linear subspace, but we do not introduce a
topology in it. We notice that the nulI subspace N p is well defined, if (10.57) is
interpreted in the sense of distributions [41]. Also the bilinear form A: D ~ D*,
given by (10.59); however, the operator P: D ~ D*, given by (10.58) is not. Many
technical difficulties are avoided by leaving the operator P out of the discussion.
TREFFIZ METHOD 237
N A= {U E D j u= :~ = 0, on OQ} (10.60)
Due to (10.60), the relevant boundary values for Laplace and reduced wave
au
equations (10.60), will be u and a;; , on oQ. We notice that given u E D and v E D,
one has that au ov
u = V' -= - on oQ (10.61)
, on on'
if and only if u - vENA; i.e. two functions u E D and v ED have the same
relevant boundary values, if and only if, u - vENA.
Similar notions can be applied to any linear differential equation. Let us
consider the biharmonic equation:
L1 2 u=0; on Q (10.62)
which occurs, for example, in connection with incompressible flows at low
Reynolds numbers. Define
<Pu,v)=SvL1 2 udx (10.63)
o
Then
oL1u ov au OL1V}
<A u, v)= S {v ---L1u-+L1v - - u - - dx (10.64)
00 on on on on
Again, a convenient definition of the space D is (see equation (10.4)):
Then, A as given by (10.64) is well defined, and N p can be taken as the linear
subspace of D which satisfies (10.62) in the sense of distributions. The operator
P: D -> D*, given (10.63), is not defined for this space D, and we leave it out from
our discussion.
The null subspace N A , is
au oL1u
NA= {uED u=a;;=L1u=Tn=O'
j
(I 0.66)
where v is the viscosity. In this case, it is convenient to define the bilinear form
P: D -> D*, by
<P 11, v) = S {v· (vL1u - V p) + qV· u} dx (10.68)
o
238 TREFFIZ METHOD
Here, u stands for a pair of functions; u which is vector valued and defined in Q,
and p scalar valued and also defined in Q. With V, we have associated the pair v, q.
Then
(10.69)
Elements of the linear space D will be pairs u= [u, p] such that the traces u and
v ~ - p n are well defined and span HO (oQ). One must also require that the set of
on
functions N pC D which satisfy Stokes equations (10.67) in the sense of distribu-
tions be well defined. In general, P: D ---+ D* may not be defined in this space. The
null subspace
{u
NA = E Diu = v :: - p n = 0, on OQ} . (10.70)
12 = {u E D I :: = 0, on OQ} (l0.76b)
and
(10.76 c)
I 2 ={UEDlu=Llu=0,onoQ} (10.77 b)
and
ou
{U E D I a;; oLl u }
13= = 8;;= 0, on oQ (10.77 c)
It has been shown [28, 34] that when {II, /z} is a canonical decomposition of D,
then II and /z are necessarily, completely regular and
(10.81)
Now, condition (10.80) is equivalent to the requirement that given any U E D, one
can find elements UI E II and U2 E I z such that
U = UI + Uz (10.82)
In the presence of equation (10.81), this representation of u, is unique except for
elements of subspace N A; more precisely, if ui E II and U2 E /z are such that
U= ui + U2 (10.83)
then UI - ui E NA and U2 - U2 E N A. Taking into account that NA is the set of func-
tions with vanishing boundary values, it is seen that the boundary values of UI and
Uz are uniquely defined. Thus, when a canonical decomposition {II, /z} is
available, representation (10.82) supplies a convenient manner of dividing the
information on the boundary values of the function U into two parts, UI E II and
Uz E /Z, which is useful in the formulation of many boundary value problems.
For the reduced wave equation, the pair {II, /z}, defined by (l0.76a) and
(10.76 b), constitutes a canonical decomposition of the space D, with respect to A,
as defined by (10.59). In this case, the representation (10.82), breaks the boundary
information in the following manner
OU OU I
U - U . - = _. on oQ (10.84)
- z, on on'
The pair {II, l]l, given by (l0.76a) and (l0.76c), is also a canonical decomposi-
'*
tion, whenever rJ. O. In this case, if U = UI + U3, with UI E II and U3 E h, then the
boundary values are given by
OU ou I OU3
U=UI+U3· - = - + _ . on oQ (10.85)
, on on on'
Ifwe define
14 = {U E Diy ~: + (ju = 0, on OQ} (10.86)
it is easy to see that {I3, I4l is a canonical decomposition, whenever rJ. (j - fJ y O. '*
Clearly, the previous ones are particular cases of this more general canonical
decomposition.
For the biharmonic equation, the following pair is a canonical decomposition
12 = {a E DI v ~: - p n = 0, on OQ} (10.89 b)
(10.97)
In this section general examples of Green's formulas are presented. Many of the
operators listed are formally symmetric in the classical sense; others can be
included due to the extension of this concept introduced in the algebraic theory of
boundary value problems [19, 34] which supplies the basic frame-work for this
chapter.
Elliptic Equations
This subject is classical. The reader is referred to the book by Lions and Magenes
[41]. The extension of such formulas to problems with prescribed jumps can be
done along the lines presented in Section 5. A general discussion of Green's
formulas from the point of view of the algebraic theory will appear soon [19].
Xz
Fig. 10.2
U * V = S u (T - t) v (t) dt (10.100)
o
is used. Let A = P - P*, then
A=B-B* (10.101)
where
ov
<B u, v) = S u * -;- dx - S v (T) u (0) dx (10.102)
oQ un Q
(ii) The Wave Equation. The incorporation of this equation in the frame-work here
presented is similar.
Taking the region Q x [0, T] and the linear space of functions D, as explained
before, define P : D ---> D* by
o2u )
<Pu,v)=Sv* ( -2-Llu dx (10.103)
Q ot
with the convention (10.100). A Green's formula for this operator is obtained
taking B : D ---> D* as
<B u, v) =
~
~v
Ju * ~ dx - J{v (T) ~
Q ~
(0) + ov (T) U CO)}
~
dx (10.104)
Formulas (10.102) and (10.104) are suitable for application to initial value
problems when the function u is prescribed on the lateral boundary of the space-
time cylinder (Fig. 10.2). More general boundary conditions can be treated by
using the Green's formula of the Laplace operator, associated with the canonical
:iecomposition defined by (10.76 c) and (10.86).
Elasticity
Let the elastic tensor C ijpq be Coo (Q), satisfy the usual symmetry conditions [48]
(10.105)
md be strongly elliptic; i.e.
Cijpq C,i nj c,p I]q > 0 whenever I c, I =l= 0, III] I =l= 0 (10.106)
244 lREFFIZ METHOD
(i) Static and Periodic Motions. Let D = HS(Q) = HS(Q) EB HS(Q) EB HS(Q), s ~ 2.
Define OU p
Tij (u) = Cijpq -;- , on Q (10.107)
uX q
(10.108)
Then, A = P - p* is given by
<A u, v) = S {Vi Ti(U) - UJi(U)} dx (10.110)
oQ
where
Ti(U) = Tij(U) nj. (10.111)
An operator B : D -> D* that decomposes A is
<Bu,v)=- S Ui Ti(V) dx. (10.112)
oQ
There are many more.
(ii) Dynamics. Let D be a suitable linear space of functions defined on Q x [0, T].
Define ( 02u. )
<Pu,v)=lvi* Q ot 2' -J'iU dx (10.113)
oV· ou· }
<Bu,v)= S ui*Ti(v)dx-S {Ui(O)~(T)--' (O)Vi(T) dx (10.115)
oQ Q ut ot
Such problems occur in many applications. In potential theory, for example, the
jumps of the function and its normal derivative are usually prescribed, while in
Elasticity the prescribed functions are the jumps of the displacements and the
tractions. Variational principles for some of these problems were developed by
Prager [49] and Nemat-Nasser [50,51] presented more recent surveys. Here, general
Green's formulas for such problems are developed which are applicable irrespec-
tively of the specific operators.
a'R 0 a'[
a"[
a"R
R [
Fig. 10.3
Consider two neighboring regions Rand E (Fig. 10.3), let 0' R = 0' E be the
common boundary separating them; in addition, 0" Rand 0" E will be the
remaining parts of the boundaries of Rand E, respectively. Let DR and DEbe two
linear spaces; in the applications to be made their elements will be functions
defined on Rand E, respectively. Consider the product space 15 = DR E8 DE;
elements U E 15 are pairs U = {UR' UE} where UR E DR while UE E DE. Given
operators P R: DR -4 D'k and PE: DE -4 D'k, define P: 15 -4 15* by
(10.166)
This additive property is usually satisfied when P: 15 -4 15* is defined by means
of an integral on the region RuE. From (10.116) it follows that
(10.117)
The symbol IVA will be used for the null subspace of 1: 15 -4 J5 -4 15*. A linear
subspace S c 15 will be considered. Elements U = {UR' UE} E S will be said to be
smooth. When U = {UR' UE} is smooth, UR E DR and UE E DE will be said to be
smooth extension of each other.
Let S c J5 = DR E8 DEbe a linear subspace. Then S is said to be a smoothness
relation if every UR E DR possesses at least one smooth extension UE E DE and
conversely. A smoolnness relation S is said to be regular or completely regular for
P, when as a subspace, it is regular or completely regular for P, respectively.
Therefore, a smoothness relation S is regular when
a) (lO.IISa)
and
b) (1 u, 0) = 0 Y UE S .and 0E S (IO.IISb)
Similarly, it is completely regular when
<1 u, 0) = <AR UR, VR) + <AE UE, VE) = 0 YO E S ¢;> UE S (10.119)
246 TREFFIZ METHOD
This is called jump operator [19,34] because it characterizes the jumps since
(10.126)
../u=V·kVu+Qu (10.127)
(PRUR, VR) = r ov OU
v../u dx + a!R uk a;;dx - aIR v k a;;dx (10.128)
(10.129)
Observe that the unit normal vector n is taken pointing outwards from the region
of integration. Equation (10.129) implies that
(10.131)
TREFFlZ ME1HOD 247
Given a = {UR, UE} E 15, let be a' = {UR, ue}, where UR and ue are smooth extensions
of UE and UR, respectively. Write
[a] = {[a]R, [a]e}; fi = {fiR, fie} (10.132)
Applying definitions (10.121) and (10.131), it is seen that
l
equation (10.129) reduces to
<AA
a, ( S
= { k [a] -ov
WR on
- v k -oa]} dx
A
on
(10.134)
2 <J a, 0) = S
G'R
{v lk oa ] - k ov [a]} dx
on on
(10.136)
is clear. Here
ov OU
<B
A
a, 0) = s uk - dx - S v k - dx (10.139)
G,(RuE) on G,(RuE) on
In a similar fashion for static and quasi-static elasticity, when the smoothness
criterium consists of continuity of displacements and tractions, one obtains for the
jump operator [28]
2 <J a, 0) = S {vdTj(u)] - [a;] Tj(v) } dx (10.140)
G'R
This yields corresponding Green's formulas.
The formulation of Greens's formulas in discontinuous fields here presented is
applicable to arbitrary formally symmetric operators which are linear. Thus, for
example, the biharmonic equation or Stoke's problem are included. Green's
formulas for two phases systems have also been derived in this manner [28].
248 TREFFIZ METHOD
Under very general conditions N pC Ip is T-complete for Ip [19, 28, 34]. For the
representation of solutions it is, however, of greater interest to have denumerable
subsets fjJ c N p which are T-complete. Examples of such systems are given in
Tables 10.1 and 10.2. It has interest to mention that for the reduced wave equations
the author has shown that a system of plane waves, which have a very simple
structure, is T-complete in any bounded and simply connected region [5].
In these tables I n (r) and H~I) (r) are Bessel and Hankel functions of the first class
[52, 53]. P'/, is the associated Legendre function, while jn and hh are the spherical
Bessel and Hankel functions [52]. We recall, in addition, that the T-complete
systems given in Tables I and 2 for Laplace equation in a bounded region are
harmonic polynomials expressed in polar and spherical coordinates. Observe that
the detailed shape of Q is arbitrary.
Laplace Equation
{I, rncos n e, rn sin n e} {Ln r- n cos n e, r- n sin n e}
Reduced Wave Equation L1 u + u = 0
{Jo (r),Jn (r) cos n O,Jn (r) sin n OJ {Hb l ) (r), H~l) (r) cos n 0, H(I) (r) sin n OJ
n= 1,2, ...
TREFFIZ METHOD 249
Laplace Equation
{r" P'f, (cos 0) eiqq>} {r- n - l P'f, (cos 0) eiqq>}
Reduced Wave Equation
{jn (r) P,{ (cos 0) eiqq>} {h~l) (r) P'/, (cos 0) eiqq>}
In each of these examples, one can give to Pfl, the structure of a Hilbert space.
Possible choices for the corresponding inner products are
ou ov 0 L1 u oA v}
(ii) S {uv+--+L1uAv+---- dx (l0.147b)
oQ on on on on
(iii) (10.147 c)
250 TREFFfZ METIIOD
With these inner products, the linear space!iJ is isomorphic to the following Hilbert
spaces:
(i) HO (iJQ) EB HO (iJQ) (10.l48a)
(ii) HO(iJQ) EB HO(iJQ) EB HO(iJQ) EB HO(iJQ) (10.148 b)
(iii) HO(iJQ) EB HO(iJQ) (10.148c)
Now, given any canonical decomposition {II, lz} it is possible to chose the
Hilbert-space structure so that the associated operator B: D --> D* (equation 10.94)
is given by
(10.149)
Thus, for example, when the inner product (10.147a) is used, equation (10.149)
yields the operator B associated with the canonical decomposition given by (10.76 a
and b). The same happens if this decomposition is replaced by (10.76a and c).
When one uses the inner product (10.147b), equation (10.149) supplies the
operator B: D --> D* associated with any canonical decomposition corresponding to
the biharmonic equation; for example, those given by equations (10.87) or (10.88).
For Stokes problem the inner product can be (10.147 c) and a possible canonical
decomposition is defined by (10.89).
For the formulation of the general boundary value problem to be considered here,
we assume there is a canonical decomposition {II, lz}' and an operator B: D --> D*
such that (10.92) is a Green's formula. Using the representation (10.82), we
formulate the problem as follows; find U E N p, such that
(10.150)
where VI is a given element of II.
Let JV p = N pi N A e!iJ =:#, be the linear space generated by the boundary values
of solutions of the homogeneous equation. Then every U E JVp can be written as
(10.151 )
where U1 Ef1 = 11 INA while U2 Ef2= lzIN A. Let JV 1 ef1 be the range of values
taken by u), in (10.151), when U ranges overJVp . Similarly, letJV2 e f2, be the
range of values taken by U2, in (10.151), when U ranges over JVp .
Given a system offunctions&iJ= {WI, W2, ••• } e N p, write
Wa= Wa1 + Wa2 (10.152)
We denote
(10.153)
Clearly, we will be able to approximate the boundary values of every solution of
(10.150), if and only if
(10.154)
Here, the bar refers to the closure of JV 1•
lREFFIZ MElHOD 251
A result similar to Theorem 10.1, holds in this more general context. Assume
I p = N p + NA is completely regular and f!iJ c N p. Then the following statements are
equivalent:
(i) pj c.,;Yp c *,spans.,;Yp in the metric of *';
(ii) /1& c.,;Yp is T-complete; and
(iii) spanf!iJ 1 = f l while spanf!iJ2 =f2 (10.155)
Therefore, when .?IJ is T-complete, it is possible to construct approximating
sequences N
UN = L. a~ W n ; N = 1,2, ... (10.156)
n=1
u~= L.
n=1
b~Wn' (10.159)
where the coefficients b~ satisfy, for every fixed N, the system of equations
N
(Wm2' UI) = L. b~ (Wml' Wnl)
n=1
(10.160)
References
I Zienkiewicz, O. c., The Finite Element Method in Engineering Science. McGraw-Hill,
New York 1977
2 Zienkiewicz, O. C., Kelly, D. W., Bettess, P., The coupling of the finte element method
and boundary solution procedures. Int. J. Num. Meth. Eng., 11, pp. 355-375, 1977
3 Brebbia, C. A., The Boundary Element Method/or Engineers. Pentech Press, London 1978
4 Brebbia, C. A., Boundary element methods. Proc. Third Int. Seminar, Irvine, Ca., 1981;
Proc. Fourth Int. Seminar, Southampton 1982, Springer-Verlag, Berlin-Heidelberg, New
York
5 Sfmchez-Sesma, F.J., Herrera, 1, Aviles, J., Boundary methods for elastic wave diffraction-
application to scattering of SH waves by surface irregularities. Bull. Seism. Soc. Am.,
72 (2), pp. 473-490, 1982
252 TREFFIZ MElHOD
6 Rektorys, K., Survey of Applicable Mathematics. ILIFFE Books Ltd, London 1969 .
7 Trefftz, E., Ein Gegenstuck zum Ritz'schen Verfahren, Proc. Second Int. Congress Appl.
Mech., Zurich 1926
8 Herrera, I., Boundary Methods for Fluids. Finite Elements for Fluids. IV, Gallagher, R. H.
(Ed.), John Wiley & Sons Limited, Chapter 19, pp. 403-432,1982
9 Mikhlin, S. G., Variational methods in mathematical physics. Pergamon Press, 1964
10 Rektorys, K., Variational methods in mathematics, science and engineering. D. Reidel
Pub., Co., 1977
11 Kupradze, V. D. et aI., Three dimensional problems of the mathematical theory of
elasticity and thermoelasticity. North Holland, 1979
12 Amerio, L., Sui ca1culo delle autosoluzioni dei problemi al contorno per Ie equazioni
lineari del secondo ordine di tipo ellitico. Rend. Acc. Lincei, 1, pp. 352-359 and 505-509,
1946
13 Fichera, G., Teoremi di completezza sulla frontiera di un dominio per taluni sistema di
funzioni. Ann. Mat. Pura e Appl, 27, pp. 1 - 28, 1948
14 Picone, M, Nouvi metodi risolutivi per i problemi d'integrazione delle equazioni lineari a
derivati parziali e nuova applicazione delle transform ate multi pIa di Laplace nel caso
delle equazioni a coefficienti constanti. Atti Acc. Sc. Torino, 76, pp. 413 - 426, 1940
15 Kupradze, V. D., On the approximate solution of problems in mathematical physics.
Russian Math. Surveys, 22, (2), pp. 58-108, 1967, (Uspehi Mat. Nauk, 22 (2), pp. 59-
107, 1967)
16 Vekua, I. N., New methods for solving elliptic equations. North Holland Pub. Co., 1967
17 Colton, D., Watzlawek, W., Complete families of solutions to the heat equation and
generalized heat equation in Rn. Jour. Differential Equations, 25 (I), pp. 96-107, 1977
18 Colton, D., The approximation of solutions to initial boundary value problems for para-
bolic equations in one space variable. Quart. Appl. Math., 34 (4), pp. 377 - 386, 1976
19 Herrera, I., Boundary Methods. An Algebraic Theory, Pitman Publishing Co., 1984
20 Sabina, F. J., Herrera, I., England, R., Theory of connectivity: Applications to scattering
of seismic waves. SH-wave motion, Proc. Second Int. Conference on Microzonation, San
Francisco, Ca., 1979
21 England, R., Sabina, F. J., Herrera, I., Scattering of SH-waves by surface cavities of
arbitrary shape using boundary methods. Physics of the Earth and Planetary Interiors, 21,
pp. 148-157, 1980
22 Herrera, I., Boundary methods in flow problems. Proc. Third International Conference on
Finite Elements in Flow Problems, Banff, Canada, 1, pp. 30-42, 1980 (Invited General
Lecture)
23 Herrera, I., General variational principles applicable to the hybrid element method. Proc.
Nat. Acad. Sci. USA, 74,pp. 2595-2597,1977
24 Herrera, I., Theory of connectivity for formally symmetric operators. Proc. Nat. Acad. Sci.
USA, 74,pp. 4722-4725,1977
25 Herrera, I., On the variational principles of mechanics. Trends in Applications of Pure
Mathematics to Mechanics, II, Zorsky, H. (Ed.), Pitman Publishing Limited, pp. 115-128,
1979
26 Herrera, I., Theory of connectivity: A systematic formulation of boundary element
methods. Applied Mathematical Modelling, 3, pp. 151-156, 1979
27 Herrera, I., Theory of connectivity: A unified approach to boundary methods. Variational
Methods in the Mechanics of Solids. Nemat-Nasser, S. (Ed.), Pergamon Press, Oxford and
New York, pp. 77 - 82, 1980
28 Herrera, I., Variational principles for problems with linear constraints. Prescribed jumps
and continuation type restrictions. Jour. Inst. Maths. Applics., 25, pp. 67 - 96, 1980
29 Herrera, I., Sabina, F. J., Connectivity as an alternative to boundary integral equations.
Construction of bases. Proc. Nat. Acad. Sci. USA, 75 (5), pp. 2059-2063,1978
30 Herrera, I., Boundary methods. A criterion for completeness. Proc. Nat. Acad. Sci. USA,
77 (8), pp. 4395-4398, (1980
31 Gourgeon, H., Herrera, I., Boundary methods. C-complete systems for the biharmonic
equation. Boundary Element Methods. Brebbia, C. A. (Ed.), Springer-Verlag, Berlin, pp.
431-441, 1981
TREFF1Z METHOD 253
32 Herrera, I., Gourgeon, H., Boundary methods. C-complete systems for Stokes problems.
Computer Methods in Applied Mechanics and Engineering, 30, pp. 225 - 241, 1982
33 Herrera, I., Boundary methods: Development of complete systems of solutions. Finite
Element Flow Analysis. Kawai, T. (Ed.), University of Tokyo Press, pp. 897 - 906, 1982
34 Herrera, K, An algebraic theory of boundary value problems. Kinam, 3 (2), pp. 161- 230,
1981
35 Herrera, I., Spence, D. A., Framework for biorthogonal Fourier series. Proc. Nat. Acad.
Sci. USA, 78 (12), pp. 7240-7244, 1981
36 A1duncin, G., Herrera, I., Solution of free boundary problems using C-complete systems.
Boundary Element Methods in Engineering. Brebbia, C. A. (Ed.), Springer-Verlag, Berlin,
Heidelberg, New York, pp. 34-42, 1982
37 A1duncin, G., Herrera, I., Contribution to free boundary problems using boundary
elements. Trefftz Approach, Communicaciones Tecnicas, IIMAS-UNAM, 1983 (Com-
puters Methods in Applied Mechanics and Engineering, submitted)
38 Liggett, J. A., Liu, P. L.-F., The boundary integral equation method applied to flow in
porous media. Allen and Unwin, 1982
39 Kikuchi, N., Oden, 1. T., Contact problems in elasticity. TICOM Report 79-8, 1979
40 Oden, 1. T., Kim, S. J., Interior penalty methods for finite element approximations of the
Signorini problem in elastostatics. Compo & Math. with Applics., 8, pp. 35 - 56, 1982
41 Lions, 1. L., Magenes, E., Non-Homogeneous Boundary Value Problems and Applications.
3 Volumes, Springer-Verlag, New York-Heidelberg-Berlin, 1972
42 Oden, 1. T., Reddy, 1. N., An introduction to the mathematical theory of finite elements.
Pure & Applied Mathematics Series, 1. Wiley & Sons, New York, London, Sydney,
Toronto, 1976
43 Bates, R. H. T., Analytic constraints on electro-magnetic field computations. IEEE Trans.
on Microwave Theory and Techniques, 23, pp. 605 -623, 1975
44 Abraham, R., Marsden, 1. E., Foundations of mechanics. The Benjamin Cummins
Publishing Co. Inc., pp. 159-187, 1978
45 Herrera, I., A general formulation of variational principles. Instituto de Ingenieria,
UNAM,~lO, 1974
46 Gurtin, M. E., Variational principles for linear initial value problems. Quart. Appl. Math.
22, pp. 252- 256, 1964
47 Herrera, I., Bielak, J., A simplified version of Gurtin's variational principles. Arch. Rat.
Mechanics and Analysis, 53 (2), pp. 131-149, 1974
48 Gurtin, M., The linear theory of elasticity. In: Encyclopedia of Physics. VI a/2, Springer-
Verlag, Berlin, pp. 1- 295, 1972
49 Prager, W., Variational principles of linear elastodynamics for discontinuous displace-
ments, strains and stresses, in recent progress in applied mechanics. In: The Folke-Adqvist
Volume. Broberg, B. Hult, 1. & Niordson, F. (Ed.), Almgvist and Wiksell, Stockholm,
pp. 463-474, 1967
50 Nemat-Nasser, S., General variational principles in non-linear and linear elasticity with
applications, in mechanics today. Nemat-Nasser, S. (Ed.), Pergamon, 1, pp. 214-261,
1972
51 Nemat-Nasser, S., On variational methods in finite and incremental elastic deformation
problems with discontinuous fields. Quart. Appl. Math., 30 (2), pp. 143-156,1972
52 Jackson, 1. D., Classical Electrodynamics. Wiley, New York, pp. 65, 69, 86, 541,1962
53 Morse, P., Feshbach, H., Methods of Theoretical Physics. McGraw-Hill, New York, 1953
SUBJECT INDEX
Boundary Element
Techniques
Theory and Applications in Engineering
1984. 284 figures. XIV, 464 pages.
ISBN 3-540-12484-5
Volume 2
B.Amadei
Volume 3
Computational Aspects of
Penetration Mechanics
Proceedings of the Army Research Office Workshop on
Computational Aspects of Penetration Mechanics held at
the Ballistic Research Laboratory at Aberdeen Proving
Ground, Maryland, 27-29 April, 1982
Editors:J. Chandra, J. E.Flaherty
1983. VII, 221 pages. ISBN 3-540-12634-1
Volume 4
W.S. Venturini
Boundary Element Method
in Geomechanics
1983. 114 figures. VIII, 246 pages. ISBN 3-540-12653-8
Volume 5
M.Manzoor
Volume 6
Springer-Verlag M. B. Allen III