Dual reciprocity BEM: local versus global approximation functions for diffusion, convection and other problems

Engineering Analysis with Boundary Elements 14 (1994) 349-356 © 1995 Elsevier Science I.imited Printed in Great Britain. All rights reserved ELSEVIER

09~~- 7997(94)00046 - 8

0955-7997/94/$07.00

Technical Note Dual reciprocity BEM: local versus global approximation functions for diffusion, convection and other problems Paul W. Partridge Departamento de Engenharia Civil, Universidade de Brasilia, 70910 Brasilia DF, Brazil (Received 20 May 1994; accepted 26 October 1994) The use of the global approximation functions (elements of Pascals' triangle, sine expansions and others) in the dual reciprocity boundary element method is compared to the better known local radial basis functions for convection, diffusion and other problems in which the volume integrals considered contain first and second derivatives of the problem variables, time derivatives and sums and products of functions, including nonlinear terms. It will be shown that whilst it is possible to obtain accurate solutions to the problems considered using the global functions, a successfulsolution to a given problem depends very much on the function chosen, as well as other factors.

Key words: Boundary elementmethod, dual reciprocity,approximation functions. INTRODUCTION

global approximation functions including elements of Pascal's triangle, trigonometric series and others is mentioned; however, the use of these functions has been unpopular due to the singular or nearly singular nature of the resulting matrices. Recently, it has been shown 5 that these matrices can be conveniently handled using a singular value decomposition algorithm. The radial basis functions are referred to as being 'local' due to the fact that they interpolate the function being approximated only in the neighborhood of the point in question, whereas the 'global' functions interpolate over the entire geometry. In Ref. 5 a series of linear problems in which the volume integral is a function of position, and two cases in which the volume integral is a nonderivative function of the problem variable are considered. It is shown for these cases that the global functions are accurate, and converge better than the local functions. Here, the use of the global functions will be considered for (i) diffusion, in which the volume integral is a function of Ou/Ot; (li) convection, in which the volume integral is a function of Ou/Ox, and (iii)a nonlinear material problem in which the volume integral is a nonlinear product of functions including second derivatives of the problem variable.

The dual reciprocity method (DRM) is a technique for taking domain integrals to the boundary in boundary element analysis and is used in cases where there is no fundamental solution to the complete problem governing equations or in cases where such a solution exists, but is cumbersome to handle, i.e. nonclosed form for instance. The method can also be used to deal with nonlinear terms, including those due to boundary conditions, The method was introduced by Nardini and Brebbia in 1982~ and a complete survey of applications was published in 1992.2 The DRM is by far the most successful technique introduced to date for dealing with domain integrals, and the only general method other than that of internal cells) When local radial basis approximation functions f -- 1 ÷ r ÷ r ~ ÷ . . . are employed, the method is very straightforward to implement.2 The convergence properties of the radial basis functions have been demonstrated.4 The method poses no limitations on the choice of approximating functions which can be used. In the original paper on the DRM,~ the possibility of using 349

P. W. Partridge

350

The computer implementation of the global functions will also be discussed,

solutions, fi, which, for the case of the Laplace operator, satisfy

v2 j THE DUAL RECIPROCITY METHOD

(6)

The method itself implies no limitation on the choice of

f. In Ref. 2, and in most work to date, the radial basis or Domain integrals in BEM analysis may arise due to a variety of effects including body forces, boundary conditions, nonlinear terms and others. In what follows, a Laplace operator will be considered for the sake of simplicity; however, the same arguments apply to any operator for which a fundamental solution is available, Consider V2u = b

(1)

subject to boundary conditions u = ~ on r~ and

q = Ou/On = ~ on /'2, where

1" 1 -]- 1"2 -~" 1" is the total boundary and n is the normal to the boundary; b will be considered to be a general fight-hand side and may be a function of position, time or the problem unknown, u, or any derivative, and may contain summations or products of terms, including nonlinear terms. A weighted residual formulation is applied to produce the basic integral equation, i.e.

Ia(V2u- b)u*df~ = Ir~ (q - (t)u'dF (2)

where fl is the solution domain, and u* is the fundamental solution such that V2u* = - A

f = 1 + r + r2 + . . . (7) are employed. These functions are shown to be straightforward to implement and results are obtained for a wide range of problems to at least engineering accuracy. In Ref. 5, the global functions: f = 1,x,y, x2,xy, y2... (8)

f = 1,sin(x),sin(y),sin(2x),sin(x)sin(y),...

amongst others are used. The first set consists of elements of Pascal's triangle, and the second set is a sine expansion. Note that no summation is implied in eqn (8). Substituting eqn (6) into eqn (5) one obtains

N+L b = ~ %.(V2t~y) y= 1 in such a way that

(9)

~+L V2u = ~ aj(V2~j)

t

- J,L,(u - O)q*dr

local functions

(3)

where A is the Dirac delta function. In the case of the Laplace operator q* = 0u*/0n and u* = 1/(2~r)ln(1/r) in two dimensions, where r is the usual boundary element distance function, Equation (2) is integrated twice by parts. After taking into account eqn (3) and grouping all the boundary terms together, one obtains for a point i called a source point:

(10)

~= ~ Multiplying by the fundamental solution and integrating over the domain gives

I

(V2u)u*dfl = ~

aj (V2~y)u*dfl

(11)

j=l

Integrating twice by parts for each source point i and taking eqn (3)into account leads to j j

¢iui + q~.ud1"N+L /

u~qd1"= t

~^J*^

~

(12) I

)

aj ctuo+ qtu]dF- u~(l~dF

j=l

ciui+Iruq*d1"-Jrqu*d1"-Iflbu*df~

(4)

where the constant c, which depends on the geometry at the point i, will be evaluated later, The domain integral term on the right-hand side will be taken to the boundary using the DRM; b is approximated as follows: N+L

b -~ E

Comparing eqn (12) with eqn (4) one can see that the domain integral has been expressed in boundary only terms; ~ is analogous to q*. In the case of the Laplace operator, ~ = Of~/#n. Diseretizing, and with the summation over the boundary elements, k, replacing the integrals, one obtains for a source node i: N

aJJ)

(5)

j=l

where N + L is the number of boundary plus internal points, a] are a series of unknown coefficients, andJ) are a set of approximating functions. The key to the method is the existence of a corresponding set of particular

N

ciui + ~ HikUk -- ~ , Gikqk = k=l

k=l

i/ N N ~ ~ j~__llOtjlCi~ljq-~'~I'Iik~kJ--k~l Gtk(tkj]

~V+L(

k= 1

=

J (13)

Dual reciprocity BEM Int[~A]

After appfication to all nodes, eqn (13) can be expressed in matrix form as N+L Hll -- Gq = E o~j(H~Ij -- G~ j ) j=

= E

m [ n l x m - ~ + 2y n + 2~

× usual meaning in BEM. Note been placed on the leading ~/N__l H/j, i ~ j. Collecting the matrices one can drop the

H u - Gq = ( H I ~ - G~)c~

(15)

(m - 2k + 2)!(n + 2k)t

for m < n

(20)

A similar equation exists for ~. It is necessary to calculate the powers m on x and n on y for a given j. A short code which will do this, storing the results in an array IP, is given below: C CALOJLATE H ~ID N FOR GIVEN 3 IN (19) AND (20) C NL1 TOTAL NUMBER OF BOUNDARY AND INTERNAL NODES

Ifb in eqn (1) is a known function, ~ is obtained explicitly from eqn (5), which, written in matrix form, becomes ¢x = F - l b

(--1)k+l

k=l

(14)

1

where H and 13 have the that the terms ci have diagonal of H using Hu = /~j and ~j vectors into summation, obtaining

351

J=0

(16)

DO i II=l, (I~1) DO 1 JJ=l, II

The problem is of size N x N and internal solutions are obtained separately. 2 If b is an unknown function of u then eqn (15) is written: Hu-

Gq = ( H I ~ - GI~)F-'b

COMPUTER IMPLEMENTATION The specific case of a global function based on the elements of Pascars triangle will be considered; other global approximation functions are treated in a similar way: f = 1, x, y, x 2, xy, y2, x3,... (18) a given term J) is expressed in Ref. 5 as being J) = xmy" in which case ~ will be given by 5 Int ["2-'~] E

m!n!xm+2kyn-~k+2 (m + 2k)!(n - 2k + 2)!

(19)

I P ( J , t)--H IP(J,2)=N IF (J. GE. (NL1)) GOTO2 CONTINUE CONTINUE

1 2

The F matrix (eqn (16)) will be singular or nearly singular if the global approximation functions are employed in nearly all cases. In order to invert F the routine SVDCMP from Ref. 6 can be employed as suggested in Ref. 5. The singular value decomposition technique is based on the fact that any matrix may be expressed as A = UWV (21) where U and V are orthogonal and W is diagonal, in such a way that A -~ = V r W - ~ U r

(22)

W -~ is easily found due to the diagonal nature of W. Before using eqn (22), however, the W matrix, stored as a vector, should be operated on, setting to zero those elements which are smaller than the accuracy of the computer in use. 6 W also permits the calculation of the conditioning number of the matrix being inverted. 6 It is not necessary to use SVD to solve the system equations (which will not be singular); the usual solvers may be employed for this process. If trigonometric of hyperbolic expansions are used for f , the computer implementation is similar. In these cases the calculation of fij and #~ is simplified as only one term is necessary for each. In the ease of a cosine expansion, for instance, if fj = cos(mx)cos(ny)

(--1)k+l

k= ~ ×

M=II-JJ

(17)

and the nodal values of b will be considered according to the ease in question (time derivatives, space derivatives, etc). The problem is now of size (N + L) x (N + L) and the solution includes both internal and boundary nodes. In the case of the global funetious, the D R M matrices are obtained in a slightly different way; the rows of O, 0 and F in eqn (17) depend on the coordinates of the field points k, but the columns are independent of the coordinates of collocation points j; each column of these matrices is obtained conside~rin~ a new term in the global expansion (eqn (8)). Thus, U, Q and F are all nonsymmetric matrices if global functions are employed. The local functions are based on rkj in such a way that each entry in the matrices depends explicitly on the coordinates of both points. 5

aj =

J=J+l N=JJ-t

(23)

then for m ~> n

ftj = -

cos(rex) cos(ny)

m 2 + ,2

(24)

P. W. Partridge

352

Similar expressions hold for sine, cosh and sinh expansi°ns'5" The code given above may be used in all cases to relate m, n and j.

APPLICATIONS

Coll~qiO[I A simple convective problem is considered in Ref. 2: V2 u

=

_

Ou

Ox

(25)

b in eqn (17) will thus be the nodal values of -Ou/Ox. This is done by writing e = F/~I and differentiating: Ou 0F

~-°x~l but since/~ = F-lu, then Ou OF ~xx = 0x F-lu

(26)

When eqn (26) is substituted into eqn (17) a matrix ¢quation: HaGq = Su (27) is obtained where S = (HfJ - G ~ ) F -t ~xFF-~

(28)

Equation (27) may easily be solved for u. Note that the columns of the OF/Ox matrix are given by

Of

~xx = 0, 1,0,2x,y,... in the case of the global Pascal's triangle expansion, The geometry used for comparing the results for local and global approximation functions is shown in Fig. 1, where the ellipse in question has a semimajor axis of length 2 and a semiminor axis of length 1. The condition u = e-x (29)

x u17 uls u19 u20 u21 u22 U23 ?/29 U30 u31 u32 u33

Table 1. llem~forV2u=-O~/Ox: L = I 7 y P.A f = l + r f = l + r 3 eqn(29)

1.5 1"2 0"6 0"0 -0"6 -1.2 --1"5 0"9 0"3 0.0 -0.3 -0.9

0-0 -0"35 -0"45 -0"45 -0"45 -0"35 0"0 0"0 0"0 0.0 0.0 0.0

0.223 0"302 0"549 1-001 1"823 3.324 4"474 0"407 0"741 1.001 1-351 2.461

0.229 0"307 0"555 1"003 1"819 3.323 4"489 0"411 0"745 1-002 1.348 2.448

0.222 0.299 0"547 1"001 1.828 3.330 4"473 0"405 0"740 1.001 1.352 2.465

0.223 0"301 0'549 1"000 1"822 3"320 4"482 0"406 0"741 1.000 1.350 2.460

is imposed on the boundary in such a way that internal results for u will also be given by eqn (29) as this is a particular solution to eqn (25). Results for the approximating functions f = 1 + r, f -- 1 ÷ r 3 and f = 1, x,y, x 2, xy, y2,.., are given in Table 1 for L = 17. The f = 1 + r solution is from Ref. 2. The global expansion solution is in the column 'P A'. It should be noted that whilst the results produced by the global Pascal's triangle expansion are more accurate than those produced by the local functions, the results for f = 1 ÷ r 3 are only slightly inferior. All solutions are well within engineering accuracy. In this case the conditioning number on the F matrix for the global expansion was of the order 10l°, being necessary to zero only two elements of W (eqn (21)) using single precision arithmetic on a 386 computer. A cosine expansion: f = i,cos(x),cos(y),cos(2x),cos(x) cos(y),... (30) was tried next. The matrix conditioning number is now 1028 and all but 12 terms of W have to be zeroed. The results produced for L = 17 are omitted as errors of up to 50% were observed. Further tests were carried out using a smaller number of internal nodes. Setting L = L~ + L2 where L~ internal collocation nodes are used with the boundary solution, results at the remaining L2 points may be obtained using N N N

u i = - E HikUk ÷ E Gikqk ÷ E HikVk k=

I

k - I

k = I

N

-- E GikWk ÷ Vi

(31)

k=l

~

where v = l~Ja, w = (~a and a = F-lb; b is defined according_~to the problem in question. In this case: b= -

Fig. 1. Geometry for conveetionproblem,

F-in

(32)

i = N + Lt + 1,N + Ll + L2; ii is the solution already obtained at points N + L~. Results for the following cases are shown in Table 2: global Pascal's triangle, Ll = 1 and L~ = 5; local

353

Dual reciprocity B E M

Table 2. R e ~ d ¢ l e t V2u = - 0 ~ / 0 x ; Ll = 1 t e d L~ = 5 x y p.A p.A f = 1 + r 3 Exact (Lt = 1) (LI = 5) (LI = 5) eqn (29) ut7 u~s u~9 U2o

u21 u22 u23 u29 u30

1.5 1.2 0.6 0.0 -0"6 -1.2 - 1.5 0.9 0.3

0.0 -0.35 -0.45 -0.45 -0"45 -0"35 0.0 0.0 0-0

0.240 0.329 0.615 1.119 1"953 3"398 4.507 0.467 0.864

0.226 0.294 0.535* 0.995 1"825" 3.324 4.474 0.391 0.729

0.224 0.301 0.550* 1.006 1"835" 3.328 4.470 0.408 0.742

0.223 0.301 0.549 1.000 1"822 3.320 4.482 0.406 0.741

u32 -0"3

0"0

1.523

1"350

1-359

1.350

u33

0"0

2-608

2"464

2"470

2.460

-0"9

one can proceed to a n e w timestep. This process can be repeated as many times as desired. The problem considered is that o f a square plate at 30°C cooled by the application o f a thermal shock (u = 0 all over the boundary). This problem was originally studied using finite elements by Brnch and Zyvoloski, 7 who reported an exact solution: oo oo - , ~ - . , _A n sm . mrx sm . mry u = ~~_~Z~, ~ Ly n = 1j= t (37)

[ {rxn2 2

]

where

f=l+r3, L ~ = 5 . T h e L~ internal nodes used with the boundary solution are marked '*'; note that due to symmetry only 12 of the 17 results are listed. With L~ = 1 the global Pascal's triangle results are not accurate. With L~ = 5 the average error in the global function solution is 0.8%; in the case o f the local solution the average error is 0"4%, thus both solutions are perfectly acceptable, however better accuracy for the global function is not observed in this case.

Difft~ion In this case the isotropic equation: V2u = k1 0u 0t

(33)

will be solved. Using D R M (eqn (17)) becomes 2 1 ^ 0u H e - G q = ~ ( H U - G ( ~ ) F -~ ~'~

(34)

The nodal values o f the time derivative are computed using a simple finite difference approximation: 0u 1 ~--[ = A-~(ut+at - ut)

(35)

where A t is the timestep. Putting C = 1 / k ( l - I ~ G O ) F -~ and introducing u = (ut+at + ut)/2 and q = qt+At the final matrix equation becomes ( 2 C ) ~ +H

(2C-H)ut ut+At- 2~lt+A t =

~'~

A, = n ~ [ ( - 1 ) ~

- 1][(-1)~ - 1]

The problem geometry is shown in Fig. 2. In all analyses, the numerical values adopted were L~ = Ly = 3, K~ = Ky = 1.25 and u0 = 30. This problem was solved using At = 0"05. For this case a global sine expansion seems appropilate, given the nature of the exact solution (eqn (37)). The results shown in Table 3 are for 33 internal nodes at time t = 1.2; the node numbers are from Fig. 2. In this case the global and local function solutions are seen to be equally accurate, and both are much better than the F E M solution, which used quadratic as opposed to linear elements. For the global sine expansion the F matrix conditioning number was 109 in such a way that only three terms o f W had to be zeroed. Single precision ailthmetic was used. Further tests, however, were run using a global Pascal's triangle expansion and a global eosh expansion, noting that whilst the latter is proposed in Ref. 5, no results are presented. In the first case the problem as posed above (N = 40, L = 33) produced a message 'No convergence in 30 iterations' from the singular value decomposition algorithm. In the case of the global cosh expansion the conditioning number on F is 102°. The attempt was completely unsuccessful, the results showing what appears to be an instability from the very first time

(36) •

which can be easily solved for Ut + A t . Setting I!t

:

=

=

:

=

:

...... •

Table 3. Diffmion problem: results at t = 1.2; L = 33

u56 u53 uT0 u69 u73

=

:

=

-

lit+At

x

y

FEM 7

f = 1 + r~

Global sine

Exact7

2.4 2.4 1.8 1"8 1"5

1.5 2.4 1.5 1"8 1"5

1-139 0.670 1.843 1"753 1"938

1.054 0-614 1.714 1-628 1-804

1.047 0.623 1.710 1"642 1.800

1.065 0.626 1.723 1"639 1-812

53' •

•

• .69.

".

.'73"7~. . 5~.

" . . . . . . .

•

"

-~. . . . . . . . . . Fig. 2. DRM discretization of diffusion problem.

354

P. W. Partridge

steps, originating at the nodes with the largest x and y coordinates.

(17) is now

b= -

U - ~ + -~--~J F

A NONLINEAR MATERIAL PROBLEM

OF~ -1 ) Ux ~xx + Uy ~yy~F ~u OF

This type of problem in heat transfer is usually solved using a Kirchhoff transformation, s This transformation eliminates the domain integrals which would otherwise arise and linearize the problem. However, if a 'third kind' boundary condition is present, even using the Kirchhoff transform the solution has to be iterative. If DRM is employed the solution is iterative irrespective of boundary condition type. The study of the nonlinear material problem was historically important in DRM as it produced a means of dealing with second derivatives of problem variables later used successfully for other problems such as anisotropic elasticity.9 A DRM solution to the case in question is given in Ref. 2 for the local exp~sion f = 1 + r 2 + r ~. The equation considered is 0 Ou 0 Ou -~x{k(u)-~xI+-~y{k(u)-~y}=O

(38)

subject to u = ~ on F 1 q=

k0U= On # on F2

(39)

q = h(uI - u) on I~3 In the above equations k is a known function of temperature, u, h is a heat transfer coefficient and uf is the temperature of a surrounding medium. Consider the case of k = k0(1 + flu) in eqn (38). In order to obtain an expression which can be treated with DRM it is necessary to isolate a Laplace operator on the left-hand side and any other expressions on the fight, There are many ways of doing this in the present case; one such result, given in Ref. 2, is

+

where U is a diagonal matrix containing nodal valves of l / f u, and Ux is a diagonal matrix containing nodal values of 1/u Ou/Ox, with a similar definition for Uy. All of the values of u used for the calculation of U, Ux and Uy are obtained from a previous iteration. Setting: ~

~ ~x + V~ ~0 ~y )' -F1 t + ( 0Ux

02F

(H-S,u =

G o~ (~nn)

0211 02F~_I Ox2 = ~--~1~ ii

(45)

The process is started with S = 0. and the Laplace equation is solved. At each subsequent iteration the matrix S is set up using the values ofu from the iteration before. A third kind of boundary condition may be introduced from noticing that Ou h 0--~= - ~ ( u f - u) Substituting eqn (46) into eqn (45) one obtains Hu - Ge + GEu = Su

(47)

24

22

" " ~ h ( ' u , Lui ~9

20

~ :

31 lra u = 300

q=0

33 (42)

35

q -- 0

X

A similar expression may be written for the second derivative with respect to y. The b vector in eqn

(46)

In eqn (47) e is a vector containing nodal values of huf/k, and E is a diagonal matrix containing nodal values of h/k. In the above, k = (1 + fu) is evaluated using values of u from the previous iteration. Equation (47) is used only on F3; eqn (45) is used on other

(41)

Since ~ = F-lu, then

(44)

then the final matrix equation becomes

26

02u

O2/_" 02F~ -1 F

S = -(HI~ - G~)F -~ [ U ~0xZ + ~ ]

A DRM approximation can be obtained for second derivatives by differentiating (eqn (26)):

Ox =b- x

(43)

•

2 4 6 8 Fig. 3. Non-linear material problem.

Dual reciprocity B E M

355

Taltle 4. Remlt~ for m for amliaear material lmflflem Node

2 4 6 8 20 22 24 26 29 31 33 35

Local

Local

Global

f = 1 + r3

f = 1 + r3

sine

(L = 25)

(L = O)

(L = 25)

306.76 305"79 304"10 301"89 307"00 311"71 314"47 315.93 313.30 310.31 308.31 307.18

309.69 308"34 305"99 302"92 307"92 313"45 316"86 318.73 316-15 313-15 311-17 310-07

308-52 306"70 304"66 302"35 306"85 311"75 314"60 316.12 313.64 310.86 309.23 308.80

~

0.tm

Kirchoff transform 306.79 305"79 304"08 301 "84 306"93 311"61 314"33 315.71 313.14 310.22 308.28 307.25

The problem studied is shown in Fig. 3. For this problem, k = k0(1 ÷/3u), k0 = 1 and /3 = 0"3. For the third kind of boundary condition h was taken to be l0 and uf---500; the remaining boundary conditions are shown in Fig. 3. The local approximation function solution used f = 1 4- r 3 and 25 evenly spaced internal nodes. The solution converged in five iterations and is given in Table 4, together with a Kirchoff transform solution, 2 which will be considered to be exact as it is obtained using a 'pure' BEM technique, without domain integral approximation. A solution using a global sine expansion is also given for the same number o f internal nodes. In addition, a local f = 1 4- r 3 result is presented for no internal nodes. The conditioning number for the F matrix for the global solution is 102~; 16 terms of W were zeroed, Serious difficulties were encountered in obtaining this solution. An over-relaxation technique had to be employed, which is not necessary for the local functions. Pascal's triangle and cosine expansions for various different numbers o f internal nodes were tried unsuccessfully. The solution shown took 21 iterations to converge, F o r L = 25 the local f = 1 + r 3 expansion is more accurate, the largest error in the local solution being 0.06%, and the largest error in the global solution being 0.5% (10 times larger). Both solutions are, however, well within engineering accuracy, It was not possible to obtain a global sine expansion solution for the case of no internal nodes; the results do not converge. The error in the l o c a l f = 1 + r ~ solution, shown in Table 4, in no case exceeds 1%.

(]en~r Point

1.4m

Y x Fig. 4. The NAFEMS benchmark problem. employed. The governing equation is 1 000 000 V2u = 52

(48)

The temperature u is maintained at 0 on the boundary. In Ref. 5, using a global Pascal's triangle expansion, a solution of 310-0 is obtained. The exact solution l° is 310.1. Further tests on the same problem showed that (i) using f = 1 4-r and one internal point at the center, 311.8 can be obtained; (ii) employing f = 1 4- r 2 4- r 3 and two internal points (one at the center and the other close by) a solution of 310"19 can be obtained. This solution is o f the same order of accuracy as that given by the global expansion cited.

CONCLUSIONS In this paper the use of the radial basis functions f = 1 4 . r 4 . r 2 4. . . . called local was compared to the use of the global functions based on elements of the Pascal's triangle, sine expansions and others as D R M approximation functions. Cases involving domain integrals containing first and second derivatives of the problem variables, time derivatives, sum and product terms, including nonlinear terms were considered. It can be concluded that: (i) For the cases where a global function solution may be found, that solution may be more accurate than that given by the radial basis functions; however, it must be remembered that it is possible to improve radial basis function solutions by choosing different f expansions. f = 1 + r 3 is frequently better t h a n f = 1 + r, for example. (ii) The existence of an acceptable global approximation function solution for a given problem depends on various factors, principally:

THE NAFEMS

BENCHMARK

PROBLEM

This problem, the geometry o f which is shown in Fig. 4, was solved in Ref. 2 u s i n g f = 1 + r obtaining a solution o f 314"6°C at the center node (15 internal nodes employed). One linear boundary element per 0-1 m is

• the function chosen • problem geometry and size The conditioning number on F does not depend on the governing equation. The acceptability of the solution seems to depend on the

356

P. W. Partridge

number of elements of W to be zeroed; if this number is proportionately large, success is unlikely. This means that a given global function may work for one geometry and not for another, for the same governing equation, (iii) For the convection problem, it was found that the local and global function solutions are equally accurate with five internal nodes, though the global functions are more accurate with 17. In the case of the diffusion problem, using 33 internal nodes, both produced equally accurate answers. Difficulties in obtaining convergence needing special techniques were encountered in the case of the nonlinear material problem considered; using 25 internal nodes the results using the radial basis function f = 1 + r 3 were superior to those obtained with the global function. It was not possible to obtain a global function solution with no internal nodes; for the same situation the local functions produced less than 1% error. Thus, on the basis of the problems considered here, the global functions are not in general superior to the local functions, though they may produce better results in a specific case. (iv) The local functions are more straightforward to program, not needing singular value decomposition; this latter cannot be expected to resolve all cases, as problems were detected with several matrices encountered. The use of the singular value decomposition algorithm depends up to a certain extent on computer architecture; 6 the larger the machine precision, the less elements of W will need to be zeroed and the larger the matrix conditioning number that can be handled, This author suggests that i r a new problem is tried using DRM, that the radial basis functions be used first. If this solution is unsatisfactory,

and here it should be remembered that all solutions as yet reported have been at least within engineering accuracy, then the global functions should be tried, and tests should be carded out to determine which is the most appropriate for a given case. The most approo priate will usually be the one that produces the smallest conditioning number o n F .

REFERENCES 1. Nardini, D. & Brebbia, C. A., A new approach for free vibration analysis using boundary elements. In Boundary Element Methods in Engineering, ed. C. A. Brebbia, Spdnger-Verlag, Berlin and New York, 1982, pp. 312-26. 2. Partridge, P. W., Brebbia, C. A. & Wrobel, L. C., The Dual Reciprocity Boundary Element Method, Computational Mechanics Publications, Southampton, Elsevier, Boston, 1992. 3. Brebbia; C. A., The Boundary Element Method for Engineers, Pentech Press, London, 1978. 4. Yamada, T., Wrobel, L. C. & Power, H., On the convergenceof the dual reciprocity boundary element method. Engng Anal. with Boundary Elements (in press). 5. Cheng, A. H-D., Lafe, O. & Grilli, S., Dual reciprocity BEM based on global interpolation functions. Engng. Anal. with Boundary Elements (in press). 6. Press, W. H., Flaunery, B. P., Teukolsky, S. A. & Vetterling, W. T., Numerical Recipes, Cambridge University Press, 1989. 7. Bruch, J.C. & Zyvolski, G., Transient two dimensional heat conduction problems solved by the finite element method. IJNME, 1974, 8(3) 481-94. 8. Bialecki, R. & Nowak, A. J., Boundary value problems for non-linear material and non-linear boundary conditions. Appl. Math. Modell., 1981, 5, 417-21. 9. Schclar,N. A. & Partridge, P. W., 3D anisotropic elasticity with BEM using the isotropic fundamental solution. Engng Anal. with Boundary Elements, 1993, 11, 137-44. 10. Cameron, A. D., Casey, J. A. & Simpson, G. B., BenchmarkTests for Thermal Analysis, NAFEMS Publications, Glasgow, 1986.