Numerical implementation of the boundary element method for two dimensional transient scalar wave propagation problems

Numerical implementation of the boundary element method for two dimensional transient scalar wave propagation problems W. J. M a n s u r * a n d C. A. B r e b b i a

Department of Civil Engineering, University of Southampton, Southampton S095NH, UK (Received May 1982)

This paper presents the numerical formulation of the boundary element method to solve two-dimensional transient problems governed by the scalar wave equation. Constant and linear time and space dependent interpolation functions are implemented and their various combinations are compared. Several examples are presented to illustrate how the time stepping technique was implemented and to show the accuracy of the solutions. Important conclusions regarding the size of time steps versus accuracy and stability can be inferred from the results. Key words: integral equation, transient numerical solution, wave propagation Introduction A comprehensive review of the boundary integral formulation used to solve transient wave propagation problems can be found in the paper by Mansur and Brebbia, ~ together with a discussion of the possible alternative numerical procedures. The present paper concentrates on the numerical implementation of the boundary integral procedures shown in Mansur and Brebbia 1 with particular reference to the twodimensional scalar wave equation. The problem is solved in a stepwise fashion and different time and space interpolation functions are used for the potential and its derivative. The initial conditions are taken into account by integrating over internal cells. The various combinations of interpolation functions are critically compared and several conclusions are presented. Numerical implementation of source density is also considered, but no application is carried out here. It is important to point out that a boundary integral equation which applies for three-dimensional analysis of the same problem has been well known for a long time 2,a and has been applied using time stepping algorithm. 4's * On leave from Federal University of Rio de Janeiro, Brazil. 0307-904X/82/040299-08/$03.00 © 1982 Butterworth & Co. (Publishers) Ltd

The paper first reviews the integral formulation of the scalar wave equation and then deals with its numerical implementation. Some problems are studied to illustrate how the time stepping technique was implemented and to show the accuracy of the numerical solution. Some important conclusions regarding the size of the time step versus accuracy and stability can be inferred from the applications. B o u n d a r y integral f o r m u l a t i o n The equation under study is the inhomogeneous scalar wave equation: V2u

1 02u - - = -- 7 c 2 /)t 2

(1)

where u is the potential, 3' is the source density, t is time and c is the velocity of wave propagation. On the P boundary the following conditions are preseribed: Essential conditions

u = a on FI

Natural conditions

au q = ~ = ?/on P2

(2)

where F = F I + P2, q is the normal derivative o f u and 0

Appl. Math. Modelling, 1982, Vol. 6, August 299

BEM for 2D transient waves: W. J. Mansur and C. A. Brebbia and ff indicate respectively, known values o f q and u on V. q can also be written as q = it • Vu, where Vu is the gradient of u, i} is the unit outward vector normal to F' and '." stands for dot product. In addition, the problem presents the following initial conditions in the domain:

The source intensity term 3', can in some cases be represented by a dirac delta function 6, such that: 7(r, t) = fs(t) 6 (r -- rs) where r s gives the position of the source. The last term in equation (4) then becomes: t+

//0 = H0

0o =

= Vo

(3)

Do and ~o are prescribed values o f u and au/Ot at t = 0. Applying weighted residuals to the above equation (for details see reference 6) one can obtain the boundary integral equation for the two-dimensional scalar wave equationt:

(8)

J". o

2c'

is(r)~/c=(t _ r) =_ R2s H[c(t -- r) -- Rs] dr

(9)

where R s = [rs-- X]. A detailed study of the problem described in this section and of how the integral equation (4) is obtained can be found in reference 1.

t+

47rChU(~k, t) =

f

J ' q u " dV dr

o

F

Numerical formulation Consider a set of discrete points r/(nodes) on the I" boundary,] = 1 , . . . , J a n d a set of values of time tn, I1 = 1 , . . . , N. u(r, t) and q(r, t) (r E I') can be approximated by using a set of interpolation functions as indicated below:

t+

+

uB+-u* •

•

0

F

dPdr

C

J

N

.0, t) = ~

1 ( VoU~"dg2

j=l

E : ' ( n ,~/(,) .~"

t?l=l

-F C 2 .

J

N

~

~

f2

q(r,t)=

1,'/ aUo u*k + e- .,I dI-uoBoa + -aR - ~ u* +,,o ~-}

rPi (ri) = 8i/

t+

ui(ri) = 8i/

(4)

4~"'(t.) = 6m,,

I2

orn(tn) = 8,, m

where: 2C

2c[c(t--r)--R] a = a(r, t/x, r) = x/[c=( t - r) = _ R ~ ] 3

H[c(t

,,~" =

--r)--R]

(S) In the expression above u* is the Green's for the infinite space and H(...) stands for the Heaviside function. In expression (4), the terms with the '0' subscript refer to initial state (r = 0), i.e.: Uo = u*(r, t/X, O)

Bo = B(r, t/X, 0) au(r, t) -

-

at

tin) (1 2)

Next, one can write equation (4) for every node 'i' and for every value of time tn, and replace the values o f u and q by their approximations given by equation (I0). The following system of algebraic equations is then obtained: N

g

y Z u,7'.7' m=lj=l N

(6)

= •

J

y. Gi';"cl}n + F}' + S n

(13)

m=l j=l

(7)

The space integrals in equation (4) must be performed with respect to the field (observation) point r and must be understood in the Cauchy principal value sense. Notice that R is the distance between r and the source point X, i.e. R = I r - - Xl, and ca is taken to be equal to 1, 0 or 1/2 for X respectively inside, outside or on a smooth part of the F' boundary.

300 Appl. Math. Modelling, 1982, Vol. 6, August

u(r/,

q}n = q(r/, tin)

and v indicates velocity as given by: -

(ll)

where 8i/is the Kronecker delta. Therefore:

u* = u*(r, t/X, r) = _ = = 2 H[c(t--r)--R] x/c ( t - - r ) - - R

v

(10)

where m a n d / r e f e r to time and space respectively, era(t), r~](r), Ore(t) and b ( r ) are chosen such that:

~2

0

O'"(t) v/(r)q~."

j=l m=l

where:

1t7/m

=

41rCi6ij6mn tn

-- f ~rT/(r) f [ cam(r)B(r,t./ri, r) P

0 i

+ 3~b'n(r.___..__)u*(r, tn/ri, r)J dr dP J ar c

(14)

BEM for 2D transient waves: 14/.J. Mansur and C. A. Brebbia tn

I"

0

l I"vou*(r, tn/ri, O) d~2

F~ = --

C2

.

£2

1

I"~uo

•

+cJ ~-u (,,t./ri,0)da a + tn

fuo B(r, t./ri, 0) dr2

(16)

R tn

S~ = i I 0

"fl~*(r tn/ri, r) d ~ dr

(l 77

~

Note that the third term on the right hand side of eqation (167 is the sum of the first and third terms of the fourth integral on the right hand side of equation (4). W]len the l" boundary is o f Kellog type 7 it can be demonstrated that ci (ci = cri) in expression (14) is given by: Ci = - -

(l 8 7

2rr

where o'i is the internal angle indicated in Figure 1. A time stepping scheme in which equation (137 is successively solved for n = 1, . . . , N can now be used to calculate unknowns u~ and q ~ at time tN. The boundary discretization adopted in this work is linear, that is, I" is represented by a series of straight line segments (elements) each one of them joining two consecutive nodes of P. In the numerical applications presented in this paper it was decided to always use (Yn(t 7 linear and try a linear or constant time interpolation function Ore(t) for the normal derivatives q. The spatial interpolation functions ~i(r7 and ui(r) were assumed constant or linear.

Q QI

\

I-

Substitution of linear or constant ¢m(t) and Ore(t) into equations (14) and (15) leads to expressions which can be integrated analytically on time. The spatial integrations were performed analytically for the case of elements with singularities and numerically for all the others, using standard Gauss quadrature formulae. In order to carry out integrations involving initial conditions (see expression (16)) the domain is discretized into triangular cells, over which initial displacements and initial velocities are linearly interpolated. The procedure used in this work is similar to that described in reference 8, where a system of polar coordinates (R, O) based at source point is used, and analytical integration with respect to the coordinate R is performed. The integrations with respect to 0 are computed numerically for each cell, using standard one-dimensional Gauss quadrature formulae. Source contribution such as shown in equation (9) can be easily considered. Notice that for them no space integration is required and that time integration is similar to the one shown in expression (15).

Applications In order to demonstrate the range of applications of the numerical procedures developed in this paper, the three following problems were studied: (i) One-dimensional rod under a Heaviside type forcing function (ii) One-dimensional rod under prescribed initial velocity and displacement (iii) Square membrane under prescribed initial velocity

(i) One-dhnensional rod under a Heaviside type Jbrchtg .function Results obtained using the two-dimensional boundary element program were compared with analytical results for a one-dimensional rod under a Heaviside type forcing function. Tire boundary element solution considered a rectangular domain with sides of length a and b (b = a/27 as indicated in Figure 2. The u displacements are assumed to be zero at x = a and their normal derivative q are also taken as null at)' = 0 and), = b for any time 't'. A t x = 0 and t = 0 a load is suddenly applied and kept constant until the end of tire analysis (E is the Young's nrodulus). Due to the topology and boundary conditions the problem is actually one-dimensional and its analytical solution can be found elsewhere. 9 Three different combinations of interpolation functions were used in the analysis as given in Table 1. The boundary was discretized into 24 constant and linear elements (Figure 3). Double nodes were used at the corners for the later model. The time was subdivided into equal intervals such that:

/3 -

C(tm -- t m - L)

= constant

(19)

/i. where l is the length of an element.

Figure 1

Internal angle used to compute c i

Combhzation 1 was tried with 13= 0.6 and gave good results for the displacments u (the order of accuracy was the same as combination 2). The numerical values of q, however, oscillated around the analytical solution, showing the onset of instability. This unstable behaviour o f q can be avoided in this particular analysis by replacing the jump

Appl. Math. Modelling, 1982, Vol. 6, August 301

BEM for 2D transient waves: W. J. Mansur and C. A. Brebbia P --H(t-O) E

q=O

P E

F u=O

q=O

/i PH ( t - O ) Figure 2

I

a = Eq

Boundary conditions and geometry definitions for one-dimensional rod under a Heaviside type forcing function

Table 1

Combination

Interpolation function

nj(r)

and

vj(r)

Linear Linear Constant

@rn(t)

Ore(t)

Linear Linear Linear

Linear Constant Constant

advantage arises from using combination 3, i.e. constant elements in space. Hence, unless otherwise state, all the BEM results presented from now on are based on combination 2. Figures 4-6 show BEM and analytical displacement results at internal and boundary points. The order of accuracy of BEM results is quite good. In Figure 7 the normal derivative of u displacements at point (a, b/2) versus ct is presented. Except for the presence of a comparatively

[

y

- - Anolybcal ........ BEM for ~3=06

E'_

Linear "11/ ( r ) a n d

v] ( r )

.i.. a

',

--

-"

',

--

I

li

=

',

-"

i

2o

3a

4a

5a

5a

7a

8a

ct

"-

Figure 4 Displacements at internal points E(a/8, b/8), F(a/2, b/2) and G(3a/4, b/2) for one-dimensional rod under a Heaviside type forcing function, rtj(r) , vj(r), ##m(t) are linear and Ore(t) is constant

I

--

=

I

Constant Figure 3

=

I

11./ ( r )

:

•

+-+

Boundary discretization for one-dimensional rod

302

Appl. Math. Modelling, 1982, Vol. 6, August

I

x

PH ( t - O ) ~l~

of the forcing function PH(t -- O) by a steep slope. It was decided, however, that this combination is not to be recommended for problems in which the u function is prescribed on parts o f the boundary. Cornbhmtions 2 and 3 were then compared to each other. It was found that for the same number of boundary elements and the same time division, better results were obtained for linear ~7/(r) and uj(r) (combination 2) than for constant ~j(r) and ui(r ) (combination 3). As the computing time is much the same for both cases it is concluded that no

a

I

and v] ( • )

~ o 1 8

~O

~

c

t

=

1

8

a

=

a a 3e a 4 2 4 Figure 5 Displacements along boundary y = 0 at times t = 0.3a/c. t = 0.9a/c, t = 1.8a/c for one-dimensional rod under a Heaviside type forcing function, rtj(r), ui(r) , dpm(t) are linear, and om(t) is c o n s t a n t ( - - ) , analytical; (e • •), BEM for 13= 0.1

BEM for 2D transient waves: W. J. Mansur and C. A. Brebbia (ii) One-dimensional rod under prescribed initial velocity

and displacement

- Anolyt~cal . . . . . . . . 8 E M for J3=06

0[ fA

For this problem the geometry and boundary conditions were identical to the previous case and, in addition, over the

B

domain f2o shown in Figure 12 the following initial conditions

I

t.la.

were prescribed:

PH ( t - O ) 20

P

a

~o

Pc

Vo(X,y) a

2a

3a

4a

ct

5a

6a

7a

=

(20)

-

E

8a

Figure 6

Displacements at b o u n d a r y points A(0, b/2), B(a/2, 0) and C(3a/4, 0) f o r one-dimensional rod under a Heaviside t y p e forcing f u n c t i o n , r/j(r), daj(r), dam(t) are linear and om(t) is constant

- -

- -

'T I

A n a l y t ~cal

Analybcol

. . . . . . . B E M for B = O 5

I

a

PH ( t - O )

. . . . . . . . B E M for [ 3 = 0 6 20

~V

v~'~l

~"

/~/~N~vl/'"

PH(t - 0 ) ~120

t,.,x . lll

J

70

8e

go

ct Figure 9 a

2a

3a

4a ct

5a

6a

7a

Normal derivative of displacement at p o i n t D(a, b / 2 ) , for one-dimensional rod under a Heaviside t y p e forcing f u n c t i o n . ~j(r), uj(r), dam(t) are linear and om(t) is constant

8a

Figure 7 Normal derivative of displacement at p o i n t D(a, b12) for one-dimensional rod under a Heaviside t y p e forcing function. Rj(r), vj(r), dam(t) are linear and om(t) is constant

a - -

Analytical

a

A n o l y t ~cal B E M for ~ = 0 8

---~)¢

I

2O

. . . . . . . .

/l J

- .......

i

PH(t-O) 20 a

20

3o

4a

5a ct

6o

70

8a

ga

100

Figure 10 Normal derivative of displacement at p o i n t D(a, b/2), for one-dimensional rod under a Heaviside t y p e forcing f u n c t i o n , rlj(r), vj(r), dam(t) are linear and om(t) is constant v-~ vrlvl/v,',ly '

a

2a

3a

4a

ct

5a

6a

7a

8a

9a

Figure 8

Normal derivative of displacement at p o i n t D(a, b/2) for one-dimensional rod under a Heaviside t y p e forcing f u n c t i o n , r~j(r), uj(r), 4fn(t) are linear and om(t) is constant

Yl

- ........

Analytical BEM f o r ~3=1 0

Px

i

small amount of noise, boundary elements and analytical solutions are in good agreement. Care must be taken on the choice of/3 in order to avoid noise, which although usually not critical for displacements, can often be excessive for tractions. In order to study the effect o f varying 13on the level of noise, four other values of /3 were studied ; 0.4, 0.5, 0.8 and 1.0 in addition to 13 = 0.6. The results for q at point (a, b/2) are plotted in Figures 8-11. One can notice that excessive noise occurred for fl < 0.6. The value 13= 0.6 was found to be the optimum for this problem.

i

o

PH ( t - 0 ) 20

2a

30

4e

5a

60

7a

8e

9a

1Oa

ct Figure 11 Normal derivative of displacement at p o i n t D(a, b/2) f o r one-dimensional rod under a Heaviside t y p e forcing f u n c t i o n , rlj(r), vj(r), rlTn(t) are linear and om(t) is constant

Appl. Math. Modelling, 1982, Vol. 6, August 303

B E M f o r 2D transient waves." W. J. Mansur and C. A. Brebbia - -

Analytical

. . . . . . . BEM for ~3=06

q=O

a

I

PH ( t - O ) 2O

/

I

a 4

q=O

Eq=PH(t -0) Figure 12 Geometry definitions, boundary and initial conditions for one-dimensional rod

a

2e

30

4e

5e

6a

7a

80

ct

\

Figure 15 Normal derivative of displacement at point D(a, b/2) for one-dimensional rod under prescribed initial conditions. ~j(r), uj(r), #fn(t) are linear and 8m(t) is constant

/I \

/ \

I I I

/ k/

'L

/

//

i \

i

\1

/

u=O

I

I

i

I

I I I Ili

i I I I I I I I I I I I

4l

\1

Figuro 13 Domain and boundary discretization for one-dimensional rod under prescribed initial conditions L,'=O

u=O

--

I

y T

F"t

/

,_.c,.~

20

Analytical

........

BEM

for

J~ =0

6

I

I I

--Point

I

I I I 111

I I I 11111

I I

u=O

.~

PH (t-O )

I

5

x

r 4

Lul~

- -

I

A

Po.n,, L'

Ak,.

I

(2

Figure 16 Geometry definition, boundary and initial conditions for membrane analysis

a

2a

3a

4a

ct

5a

6a

7a

8a

Figure 14 Displacements at boundary point A(0, b/2) and internal points 1(3a/16, b/2), H(3a]4, b/2) for one-dimensional rod under prescribed initial conditions. ~j(r), uj(r), c~m(t) are linear and om(t) is constant

The analytical solution for this problem is the same as for the previous one but with the time t dephased by a/4c, i.e.:

°) u(i)(r, t)-- lt(ii) (r, t - -~c

r~ IX I \ \/ I /\

q(n(r, t) = q0;)

V_ _ _ _ _ _ N

I/

, t-

(21)

Twenty-four linear elements were used to discretize the boundary and ~ 0 was subdivided into four triangular cells as shown in Figure 13. The time steps were such that /3 = 0.6. Displacements at points (0, b/2), (3a/16, b/2), (3a/4, b/2) and traction at point (a, b/2) are presented in Figures 14 and 15. The accuracy o f the results is similar to the one presented by the results shown in the previous problem.

304

A p p l . M a t h . M o d e l l i n g , 1 9 8 2 , V o l . 6, A u g u s t

/

71 I I

\1

Figure 17 Membrane discretized into 32 elements and four cells

BEM for 2D transient waves: Wo J. Mansur and C. A. Brebbia

(iii) Square membrane under prescribed hlitial velocity The subject of this investigation is a square membrane with an initial velocity Vo = c prescribed over the domain f2o shown in Figure 16 and zero displacements prescribed all over the boundary. The boundary was discretized into 32 elements and ~ o was divided into four cells (see Figure 17). Analytical 1° and BEM results for displacements at point (a/2, a/2) and the normal derivative of displacements at point (a, a/2) were compared. The values of u and q for/3 = 0.6 are plotted in Figures 18 and 19 respectively. Although the agreement for displacements is reasonable, it was found that a more refined time division was needed to represent q more accurately. Another Boundary Element analysis was carried out with/3 = 0.2 and the results obtained for q are plotted in Figure 20, and they show a better agreement. Finally, a last analysis was performed for which the boundary was discretized into 64 rather than 32 elements, and the value of/3 was taken as 0.6. The results for this analysis (see

o o° °iI@ t Yl

--Analytical

ol oooo ,o 0o

020 -,~

O'-

-0 20 -040 -O 60 1

I

02a

06a

I

I

a

I

I

ct

140

006 ~]

t

I

22a

Analytical ooooBEM for ~3:O6

B

O 04 -

O06L /

ool/r/

ooooo BEN1 Analyt,cal for ~=O6

Yl

00~ 0

0 02~~1~

I

18 a

Figure 20 Normal derivative of displacement at point B(a, a/2). 32 boundary elements

YT [

I

-002

-

0 ~

-

- -

-0 0 4 -O 06

-0 02V -0 04 u

O2a

- 0 06 V

I

O2a

O6 a

a

I

I

14a

I

I

1 8a

I

06a

a

14a

1

8a

2 20

ct Figure 21 Normal derivative of displacement at point B(a, a/2). 64 boundary elements

22a

ct Figure 18 Displacement at point A(a/2, a/2}. 32 boundary elements

Figure 21) were only slightly better than those for the Yl

Analytpcal

O 6 0 " ~

ooooBEM for 13=O6

previous case, apparently because unlike the rod analysis, /3 < 0.6 did not introduce too much noise into the numerical results. Conclusions

%

-0 20

60/

I O2a

I I O 0a

I a

I

1 1 14a

I 1 80

I

I 2 2a

ct Figure 19 Normal derivative of displacement at point B(a, a/2). 32 boundary elements

This paper has presented results for the analysis o f twodimensional transient scalar wave propagation problems using the boundary element method. Interpolation in time and space were used and certain recommendations regarding their combination were presented in the examples. Care should be taken in the analysis on the choice o f time intervals and boundary discretization in order to avoid contradicting the causality property 1 too far, that is, in each time step waves should not be allowed to travel between nodes far from each other. The examples presented demonstrate the remarkably good accuracy of the results and validate the application of the boundary element method to solve transient twodimensional problems governed by the scalar wave equation.

Appl. Math. Modelling, 1982, Vol. 6, August 305

BEM for 2D transient waves: W. J. Mansur and C. A. Brebbia

References 1 2 3 4 5

Mansur, W. J. and Brebbia, C. A. Appl. Math. Modelling 1982, 6, 307 Kirchhoff, G. Ann. Ph),sik. 1883.18, 663 Kirchhoff, G. 'Vorlesung Uber Math. Physik', Math. Optik, Bd. 2 (5th edn), Teubner, Leipzig, 1891 Mitzner, R. M. J. Acoustical Soc. Amer. 1967, 42, 391 Groenenboom, P. H. L. 'The application of boundary elements to steady and unsteady potential flow problems in two- and three-dimensional problems', in 'Boundary Element Methods'

306 Appl. Math. Modelling, 1982, Vol. 6, August

6 7 8 9 10

(ed. Brebbia, C. A.), Springer-Verlag, Berlin, 1981 Brebbia, C. A. and Walker, S. 'Boundary Element techniques in engineering', Butterworths, London, 1979 Kellog, O. D. 'Foundations of potential theory', Dover, New York, 1953 Telles, J. C. R. and Brebbia, C. A. AppL Math. Modellhlg 1981, 5,275 Miles, J. W. 'Modern mathematics for the engineer' (ed. Beckenbach, E. F. ), McGraw-Hill, London, 1961, pp. 82-84 Morse, P. M. and Ingard, K. V. 'Theoretical acoustics', McGrawHill, London, 1968