Boundary element method and wave equation

RESEARCH Boundary equation

NOTES

element method and wave

Djuro M. Misijenovic Mathematical Institute, Beograd, (Received Jarmary 1982)

Yugoslavia

In this paper the boundary integral expression for a one-dimensional wave equation with homogeneous boundary conditions is developed. This is done using the time dependent fundamental solution of the corresponding hyperbolic partial differential equation. The integral expression developed is a generalized function with the same form as the well-known D’Alembert formula. The derivatives of the solution and some useful invariants on the characteristics of the partial differential equation are also calculated. The boundary element method is applied to find the numerical solution. The results show excellent agreement with analytical solutions. A multi-step procedure for large time steps which can be used in the boundary element method is also described. In addition, the way in which boundary conditions are introduced during the time dependent process is explained in detail. In the Appendix the main properties of Dirac’s delta function and the Heaviside unit step function are described. Key

words:

boundary

integral

element

method,

wave

equation

Introduction The Boundary Integral Element method has become an increasingly popular method for solving problems described with partial differential equations.’ The first analytical results for boundary integral methods were developed at the beginning of this century. The method was re-examined 10 years ago using the new technology which can be efficiently applied for approximate solutions.2-4 Currently, many researchers are interested in the boundary element solution of time dependent problems. Boundary element results have been represented for parabolic partial differential equationF7 but for hyperbolic partial differential equations the theory and application are not as well developed.*-”

Formulation

of the problem

to sblve the following one-dimensional hyperbolic partial differential equation:

We want

a% ---2 at*

a5 - = 0 ax*

xE52,tER

(1)

This is the well-known wave equation which governs many problems in continuum mechanics, such as acoustical,

electromagnetic, etc. The equation presents two types of boundary conditions: ll=ii

on K2;lxR

(24

au an’4

on KJxR

(2b)

We assume that functions ii and Q are equal to zero, i.e. the boundary conditions are homogeneous. Also, equation (1) requires the following initial conditions: u = @I all at-

-

zr2

onax{

@a)

onfix{

(3b)

Here, without loss of generality, we assume the time interval as [0, T]. Ln a more general case we may have to perform shifting (translation in time domain). Also, we suppose that $2 is finite set, n = [0, L].

Developing

BIEM expression

We interpret equation (1) in terms of potential theory, i.e. u is a potential and au/an a flux along the outward-directed normal at the boundary point under consideration.

0307-904X/82/030205-04/$03.00

0 1982

Butterworth

& Co. (Publishers)

Ltd.

Appl.

Math.

Modelling,

1982,

Vol

6, June

205

Research

Note

It is easy to see that the first integral term in equation (lo), after applying (4) and the selective property of the delta function, gives:

In potential theory, the adjoint equation to equation (1) is:

a3.d - - $ a’” = -6(x at2 a2

-J,) &(r - 7)

(4)

all

11(x, f) =

which expresses the effect of unit source at the point 01, r) in the absence of other source points in the domain of interest. The asterisk (*) denotes the unbodnded character of this solution. c is the speed of propagation of the disturbance in the medium under consideration. 6(a) is described in the Appendix. It is well known that the so-called fundamental solution for equation (4) is given by:

u*=--fqct-r]

at

=+E )I

o

dy

(11)

Now, introducing the causality property of the fundamental solution, i.e.: ~(*0,,7;x,t)=O

wheneverclr-tl
(12)

and the explicit formulation of the function cl*, as given by (9) we obtain: L

t)=

-s 0

I a

--H[cr-lly--lluOt,O)dy

2~ at

L

a4.h 0) ardy

+

(6)

w dy dr = 0

at

0

u(x,

H(o) is described in the Appendix, and r is the distance between source point and point under consideration. We transform equation (1) to a corresponding integral expression using the weighted residual method, i.e. multiplying equation (1) by a function w, which gives:

all*

l.P--U--

S(

(5)

2c

0

L

(13)

Using expression (A5) and the ‘cut off’ property of the Heaviside function we have the following formulation:

0

where w is an arbitrary function to be defined. After partial differentiation of expression (6) and introducing homogeneous boundary conditions, we obtain:

11(x, t) = f [ll(x - ct, 0) + u(x + ct, O)]

(14)

As is usual in the weighting residual method, the function w distributes errors over the whole domain. The best choice for the weighting function w is to take: w=u*

u(x, T) = +(x

-cT,O)

+ u(x + cT,O)]

(8)

Function w is a function of four variables: the first two variables are the space-time coordinates of the source point and the second pair of variables correspond to the coordinates of the point under consideration. This function is given by: w=w(y,T;X,t)=-~H[c,T-t/-ly-xll

(15)

(9)

With this choice of weighting function expression (7) gives: T+E

This expression is understood in the sense of a Cauchy principal value and it allows for the solution of a much wider class of problem than the ordinary expression known as the D’A1amber-t formula. It is clear from the above that expression (14) is also valid on the upper boundary t = T of the domain of interest. There, we have:

L

a2u* at2

2az”* u dy dr ax2

--C

L

+ I( 0

au u*--u at

Math.

Modelling,

au* T+E dy = 0 at )I o

1982,

Vol

Calculation of derivatives useful invarian ts

of the solution

and some

Expression (14) (or equivalently (15)) can be used for calculating the partial derivatives of the solution with respect to time and space. They are:

(10)

Next, we integrate in time in the interval [0, T + E], E > 0.

206 Appl.

Note that all integrals are now given on the boundary (some of them are given as functional values and the rest of them disappear because of homogeneous boundary conditions).

6, June

a+ -=at

t)

c - au+ - ct, 0) + a+ 2 [ ax +-

+ ct, 0) ax

I au(x - ct, 0) + al+ + ct, 0) 2 at at

1 1

(16)

Research Note and &1(x, -=ax

t)

I au(x - ct, 0) + ~U(X + ct, 0) 2 ax ax + L

1

_ au(x - ct, 0) + au(x + ct, 0) (17) 2c [ at at

1

Note the simplicity of equations (16) and (17) which depends only on functional values of the corresponding partial derivatives. From the above expressions we can conclude that the following two quantities are constant in the direction of the characteristics:

au(x, t) c-++ ax

aLc(x, t> at

aL[(x + d, 0) + au(x + ct, 0) =c ax at au(x, t) c--ax =c

numerical results oscillate around the exact solution and the large number of points are needed to reduce this ‘noise’. The boundary element solution on the other hand did not have this problem and results were coincident with the exact ones.

Multis tep procedure A problem of practical importance is to calculate a solution with time steps as large as possible. In this case we can use the following multistep procedure. Let CT = s, where s = length of one element. We want to calculate the solution at the time t = 2T. Shifting time by one time step in expression (15) gives: u(x, 2T) = f [u(x -s, T) + u(x + s, T)]

(18) T x;s au(y, T)

+2s J ~

aflcx, t)

at

x-s

at

aucx - ct, 0) au(x - ct, 0) ax at

(19)

Equation (18) and (19) gives the characteristic directions of propagation of disturbance in the space-time domain.

dv

(21)

After introducing (1.5) (16) and (17) we find: u(x, 2T) = ; [u(x - 2s, 0) + u(x + 2s, 0)] x+2&! auk -

+f

Boundary

element method

formulation

4s

We use expressions (1 S), (16) and (17) for numerical calculations. !2 is divided in N equal parts. We describe each part with two nodes (end points) and three values-degrees of freedom (u, au/ax, au/at) at each node. Hence, we have cubic variation of the function cl in the space domain. All nodes (except the first and last) were double nodes. This approach allows the discontinuity of the space derivative of the solution to be taken into consideration. In the case for a shape that is continuous in the space domain, any discontinuity of the derivative will indicate a numerical error.

Numerical

We have solved two casesof a vibrating string with fixed ends. In both cases, the string had unit length and the speed c was chosen to be equal to 10. In the first example, the displacement is 0.1 in the middle of the string and varies linearly towards the fixed ends. In the second example the 0.1 displacement was made to act at a point located one-third of the way along the string length. These were the two initial positions of the string. Initial velocities were assumed in both cases to be zero. Numerical results with only six elements give the same values as the exact solution. The time step was given by the Courant number, K: c. At KE ---cl n ._

(20)

This time step corresponds to the time for which the disturbance advances across one whole element. If the problem is solved using finite differences the solution progress without dumping or amplification but

cm

dy

at

This formula can, by induction, be extended to: u(x, IIT) = ; [u(x -

.Y+fl.T adi 0) +-T ~ s dy, 2ns

Introducing

examples

s x-2s

0)

boundary

IIS,

0)

+

u(x

at

+

IZS, O)]

nEN

(23)

conditions

If we want to have the same formula for calculating the position of all the points in the string including the two boundary points, we need to define some values outside the string domain. In wishing to do so, we need to assume that the string extends outside R and that formula (15) can be applied in the following way. At the initial moment the integrand of (15) is equal to zero and the displacement on the left-hand side (fwed end) is equal to zero. As we have no restrictions on time in this formula, it follows for this end of the string that: u(-z,

t) = - u(z, t)

z E [0, s], t E [0, T]

(24)

In the same way, if the string passesthrough the equilibrium position u(x,O) = 0, x E [0, L] and has some nonzero velocity we can conclude for fixed end x = 0 that:

a+-z, t) auk t) =-zE[O,s],tE[O,T] at at

_

(25)

Similarconclusions are valid for the other fixed end at x=L.

Appl.

Math. Modelling,

1982, Vol 6, June

207

Research Note

Concluding

9

remarks

Brebbia, C. A. (ed.). ‘Boundary element methods’, SpringerVerIag, Heidelberg, 1981 Brebbia, C. A. (Private communication)

The problem which is under current investigation is the solution of shock waves. This is a more difficult problem as it corresponds to nonhomogeneous boundary conditions but encouraging results have already been obtained.

10

Acknowledgements

Here, we briefly describe the properties of the Dirac delta function and Heaviside unit step function which are used in this paper. These functions are of great importance in developing BIEM formulae. 6(or) is the well-known Dirac delta function with two main properties:

,

The author is grateful to Dr C. A. Brebbia, University of Southampton, for his helpful comments and support in all stages of developing these results. He also thanks the Computational Mechanics Centre, Southampton, for use of their computing facilities on VAX 1 l/750, and the Mathematical Institute, Beograd, Yugoslavia, for every support.

Appendix

6(a) = 0

a#0

(Al)

and +-

References

208 Appl. Math. Modelling,

s

6((u) do = 1

Brebbia, C. A. ‘The boundary element method for engineers’, Pentech Press, London, 1978 Brebbia, C. A. and Wrobel, L. C. ‘Recent advances in numerical methods in fluids’, VoL 1 (ed. Taylor, C. and Morgan, K.), Pineridge Press, UK, 1980, pp 1-25. Wrobel, L. C. and Brebbia, C. A. ‘Boundary elements in water resources’, NATO Instituto at Viierio. Portugal. Aunust 1981 Brebbia, d. A. and Wrobel, L. ‘Apphcations of boundary ele. ments in water resources’, 3rd Int. ConJ Fin&e Ekrnerm in Water Resour., 19-23 May, 1980, University of Mississippi, Oxford, USA, pp 2.3-2.14 Shaw, R. P. ‘An integral approach to diffusion’, Inf. J. Hear Mass Transfer 1914, 17, 693 Liggett, J. A. and Liu, P. I. F. ‘Unsteady flow in confined aquifers: a comparison of two boundary integral methods’, Water Res. Res. 1979, 15 (4), 861 Brebbia, C. A. and Walker, S. ‘Boundary element techniques in engineering’, Newnes-Butterworth, 1380 Brebbia, C. A. (ed.). Progress in boundary elements’, VOL 1, Pentech Press, London, 1981

1982, Vol 6, June

(A21

The selective property of Dirac’s function gives: +(A31

H(a) is the Heaviside unit step function defined by: H(a)=

:,

(Y>O

a<0 t and with the useful property that:

$H(o!)= 6(a)

C.44)

(A9