An efficient boundary element method for two-dimensional transient wave propagation problems

An efficient boundary element method for twodimensional transient wave propagation problems Vedat Demirel and Shen Wang Division of Applied Marine Physics, Rosenstiel School of Marine and Atmospheric University of Miumi, Miami, Florida 33149, USA (Received

Science,

October 1986; revised February 1987)

The boundary element method (BEM) has been recognized by its unique feature of requiring neither internal cells nor their associated domain integrals in the computation. The method preserves its elegance for transient problems by means of a certain time-stepping scheme that initiates the time integration always from the initial time. Unfortunately, this time-marching scheme becomes rather difficult to apply, because the computation time and storage requirement grow dramatically with the increasing number of time steps. This paper shows that a reduction of one half of the computation time as well as the storage requirement can be achieved by an efficient truncation scheme for two-dimensional transient wave propagation problems. In particular, a guiding parameter for the determination of the truncation limit is proposed, and the overall measure of the error with respect to the truncation guide parameter is established. Keywords: propagation

integral

Introduction The boundary element method (BEM) for solving steadystate problems has been appreciated for a long time because of its efficient and convenient numerical scheme. However, its extension to transient problems is relatively new, especially in the field of two-dimensional wave propagation. With BEM formulation a direct timedomain solution of the two-dimensional scalar wave equation was given by Mansur and Brebbia.’ In a subsequent paper,’ they discussed its numerical implementation and later3 gave an extension of the formulation to elastodynamical problems. As they indicated, the results are strikingly accurate, but the time-stepping scheme becomes cumbersome as the elapsing time increases. Recently, Demire14 considered a transient harbor oscillation problem and improved the technique of time stepping by introducing a truncation limit. 0307-904X/87/060411-06/$03.00 0 1987 Butterworth Publishers

equation;

efficient

transient

solution;

wave

This paper summarizes this new concept and demonstrates its advantages. With the aid of this concept, the time-marching scheme used by Mansur and Brebbia has been improved, and remarkable reductions in computation time and storage requirement have been accomplished.

Boundary element formulation The differential equation under consideration dimensional, linear, long-wave equation

is the two-

$L,2($+$$)

(1)

Here c represents the velocity of the wave propagation, and 4 is the principal unknown, say velocity potential or surface elevation, which is a function of the two horizontal

Appl.

Math. Modelling,

1987, Vol. 11, December

411

Two-dimensional

transient wave propagation problems: V. Demirel and S. Wang

dimensions x and y as well as the time variable t. Following the conventional BEM procedure, we seek an integral equation corresponding to (1). Details of this derivation are in Mansur and Brebbia.’ The resulting integral equation for wave propagation in a domain R of boundary r has the form

+

q4* dl- dz

(2) where R = 1I - [) = distance between observation

and 8 are time interpolation functions. Substituting these interpolation functions into equation (2) will yield a set of algebraic equations. One may have a different equation system for each of the two time-marching schemes, known as procedure I and procedure II. Procedure I, which is more traditional, treats each time step as a new problem. Namely, at the end of each time step the values of 4 are computed at a sufficient number of internal points in order to use them as pseudo-initial values for the next time step. In procedure II, however, the time integration process always starts from the initial time. Thus, despite the increasing number of intermediate steps as time progresses, values of 4 at internal points need not be computed. Procedure II is usually preferred because domain integrals are avoided. Substituting the interpolation functions (4) and (5) into (2) and repeating the procedure for all boundary nodes yields

point r

and source point [ The elements of the matrices [HFf], [G”“], and {OF} are

q = n. V$ = normal derivative of 4

hcf = 47&,6,,@

v = &#@t = time derivative of 4 The parameter 1 is taken equal to 1, 0, or 3 for [ inside, outside, or on a smooth part of the r boundary, respectively, and

2c = (cz(t _ # _

R2)1/2

HCc(t

-

‘)

-

R1

+q5,$dR,

2c[c(t - z) - R] (C2(t

_

92

_

~2)3/2

HCc(t

-

‘)

-

H[c(t - z) - R] = c(t - z) > R

Numerical formulation Consider an arbitrary boundary r surrounding the domain Q. Discretize the boundary into N straight segments, the domain n into M cells (say triangular subregions Q,), and subdivide the time dimension into F time steps. Furthermore, assume that functions 4 and q in equation (2) vary within each element and time step according to the space and time interpolation functions such that

drj,

tf) =

C'Y4

tf)=C%

(4)

Appl. Math. Modelling,

As emphasized above, the preferred time-stepping scheme, procedure II, requires all the time integrations to be evaluated from the initial time step to the final one. Therefore, computation time increases dramatically as the number of time steps increases. Our objective is to truncate the integration while keeping the accuracy within acceptable limits. For a given observation point r and time t,, the variation of 4* with respect to z is sketched in Figure 1. Suppose that the time integration of $* is required between the limits z = 0 and z = t,. For sufficiently large t,, there may be a moment (say tL) before which the contributions to the integral are insignificant. In other words, the integration may be started from z = t, instead of z = 0. The price for the truncation of the time integration is the introduction of an error, the magnitude of which depends on the value of t,. It follows that a criterion for setting the truncation limit t, is necessary. Now, define a parameter 6 by a simple ratio : ;i=sb”

q5* dT/sg”‘b*

dr

(10)

(5)

Here c and p are space interpolation functions, whereas y

412

(9)

Concept of truncation limit

a unit step function

The function @ is Green’s function for the twodimension wave operator. It may be considered as the effect of a source presented by an impulse at t = z located at r = [. The subscript 0 indicates initial time, z = 0. Integrals over r and over R are boundary integrals and domain integrals, respectively. All are Cauchy principal value integrals.

4(rjp

1

In equation (9), Rid represents the distance between the source point and the observation point within the subregion Qd.

with c(t - z) < R

td )

R1

1987, Vol. 11, December

As shown in (lo), the parameter

6 is not a true error

Two-dimensional

I

I

0 Figure

r=t-RR/c

s=tL 1

Variation

of qb* with

respect

transient

wave propagation

and S. Wang

c

z

to T

measure of the problem, but provides only a guide of the truncation by comparing the truncated integration with the untruncated integration of a typical function of the problem (in this case, Green’s function 4*). On the other hand, the true error measure can be described with the guide parameter 6, as will be shown later. Inserting c$* from (3) into (lo), replacing R by a constant Ravg (the average distance between the boundary segments), and defining tavgas R,,$c, we obtain the following relationship after evaluating these integrals:

Figure

Variation

of

1 :=I

tL/tF with

reipect

.

f=21f=3

to d

.

.

.

.

f=f,

.

.

.

D-t”1

p=1

p=

~_[(1_~)l(!$yl”2

[Hz1

2

8==3

[H3’] .

=l_{l+[l_(~~]1’2jl-d(~)”

V. Demirel

problems:

j’[H22] [H3;j .

[H33] . ‘\

(11)

The choice of Ravg to replace R in the derivation of (11) is optional. One may choose any other fixed distance such as &in or R,,, (minimum or maximum distance between the boundary segments, respectively). Evidently, a choice of Rmin would be relatively less accurate, whereas R,,, would not provide a large saving in the computation. Although this point deserves more attention, Ravg seems to be an acceptable parameter on all grounds. The solution of equation (11) gives a reasonable guide to the truncation limit t, for a prescribed 6. Even without rigorous error analysis, one may expect that the true overall error should be less than the guide parameter 6, since the influence coefficients of 4’s in equation (2) are proportional to the function B, which decays much faster than +*. In accordance with equation (1 l), the variation of the ratio t,/t, with respect to 6 is illustrated in Figure 2 for four different values of t,/t,,,. The fact that the ratio of t,/t, becomes higher with increasing 6 indicates that the truncation time t, is directly related to the error guide parameter 6, a large truncation resulting from higher 6. Likewise, for a given 6, t, increases with increasingly longer final time t,. Up to now, we have discussed the impact of truncation on the reduction of the computation time. Another equally important effect of this truncation is on the storage requirement. Take the influence coefficient matrices [Hpf], for instance. A list of the required matrices for each time step fl are depicted in Figure 3. As Pina’ noted, the matrices on each horizontal line are different, whereas the matrices along any diagonal line are identical. Therefore, among all the matrices in Figure 3, only the firstcolumn matrices need to be computed and stored. What

f=l

f=2

f=3

.

f=F

t truncation Figure

3

Sketch

of matrix

arrangment

the truncation contributes to this figure is shown by the shaded zone of Figure 3, where all the elements of the matrices are small. Thus, there is no need to compute nor store them. Making use of the truncation limit concept, we write equation (6) in the form

,$LIH”1{4’;

= f$

CW@f> + {@I

(12)

L

Notice that the lower limits of the summations (i.e.,f= 1) have been replaced byf, = tJAt. In other words, the lower limits (r = 0) of the time integrals in (2) have been replaced byr=t,. Let us compare the new form (12) with (6) by referring to Figure 3. Originally, a total of F(F + 1)/2 operations had to be performed under each summation sign on both sides of (6) to obtain the time history of oscillations from the initial time to the final time t,. In the new approach, however, the total number of operations is reduced to the summation of those matrices lying between the two dashed lines. In other words, the total number of operations under each summation sign is decreased by fLV;. - 1)/2. Actually, the saving in the computation time arises exactly from the reduction of the summation operations.

Appl. Math. Modelling,

1987, Vol. 11, December

413

Two- dimensional transient wave propagation problems: V. Demirel and S. Wang

Numerical examples The transient response of a rectangular harbor due to an exponentially decaying cosine wave and a train of periodic sine waves is used to display the valuable improvement of the new formulation. The domain of interest consists of two parts, domain I (open sea) and domain II (harbor); see Figure 4. The formulation and solution are conducted at each domain separately. In domain I, separating the potential into its components is conventional in linearized water-wave problems: 41 = $+ + 4, + +,, where the components &, $J,, and 4, are the incident, reflected, and scattered wave potentials, respectively. The boundary and initial conditions are specified as follows:

on the coastline

q,,=O

on

untruncated 6=0.20

Initial conditions

Boundary conditions qs = 0

0

I

DEA

I 0.0

I

2.0

4.0

6.0

8.0

’

t/T

4i + 4, + 4s = &I 4i+4r+43=-411

on DA

4,, = !$

= 0

I

For more details of the problem formulation and its analytical background, see Refs. 4 and 6. The results of the two numerical examples are shown in the following figures. The observation point was taken to be a point located at the back boundary of a rectangular harbor, which has a width-to-length ratio of 0.5 and an average arrival time tavg= 0.54. The boundary of the harbor was divided into N = 32 straight segments. Furthermore, the space interpolation functions &,p and the time interpolation function 0 were taken as constant, whereas y was taken to be linear. In these figures the value of the surface elevation, q, was nondimensionalized by (~i)max, the maximum elevation of the incident wave. Likewise, the elapsed time, t was nondimensionalized by the period of the incident wave, T Figure 5 shows the transient response due to an exponentially decaying cosine wave. The truncated results correspond to the case with guide parameter 6 = 0.20. We

Figure results

5 Comparison of truncated for an exponentially decaying

(6 = 0.20) and untruncated cosine wave input

~~ untruncated I’

6=0.20

4 (%)max

9

0

i 0

I

70.0 Figure results

6 Comparison for a continuous

I

2.0

f/T

4.0

of truncated (6 = 0.20) train of sinusoidal wave

E and untruncated input

Y

t

Figure

414

4

Geometric

Appl.

Math.

configuration

Modelling,

of a rectangular

1987,

harbor

Vol. 11, December

see that except for the presence of very small differences in the later time steps, the two curves in Figure 5 are in good agreement. Especially from the engineering point of view, these discrepancies are insignificant, since the vitally important part of the prediction, the maximum amplifications, agree almost exactly. Figure 6 shows a similar computation for a continuous train of periodic sine-wave input. As represented by the solid curve, the amplification at the back boundary increases with time until a steady state is established nearly four cycles later. As before, the truncated results with 6 = 0.20 show good agreement with that obtained by the normal untruncated procedure. On the other hand, a remarkable reduction in CPU time and storage requirement was achieved, as demonstrated in Figures 7 and 8. Figures 7 and 8 summarize required CPU time and the saving of storage capacity as a function of the number of time steps and truncation guide parameter 6. For instance, nearly a 50% saving in CPU time and a 60% saving in storage requirement are obtained for 6 = 0.20.

Two-dimensional

0

/I

+

transient

problems:

V. Demirel

and S. Wang

(Er),,,. Thus we analyzed a number of runs for the problem in Figure 6 for a fixed t, corresponding to 125 time steps but different values of 6. The variation of (Er),,, with6 is shown in Figure 9, where an approximately linear correlation between 6 and (Er),,, is exhibited. In general, for 0.10<6<0.35. The we found that (ET),,, <6/2 computations were conducted on a VAX 11/785 system using single-precision FORTRAN. As Figure 7 shows, the CPU time dramatically decreases with increasing 6. On the other hand, a choice of higher 6 implies less accuracy, as illustrated in Figure 9. In other words, the CPU time is a monotonically decreasincreasing function ing, whereas (ET),,,,, is a montonically of 6, but neither of them individually governs an optimal 6. However, we could establish an optimal choice of 6 (6,,,) by combining the contrary behaviour of the 6 dependence. For instance, we could obtain dopr for the numerical example given here by means of the sketch in Figure 10. The curve represented by triangles corresponds to the CPU values read from Figure 7 for 125 time steps, whereas the other curve is identical to Figure 9. The intersection of the two curves yields an optimal 6.

untruncated 6=0.15 6=0.20

time step Comparison of CPU times for three different choices guide parameter (5 = 0.0, 6 = 0.15, ii = 0.20

Figure 7 truncation

wave propagation

of the

/ 0

d

n

LO.15

+

6=0.20

L+ 0.00

Figure

I I

I 0

50

I

100

I

9

Correlation

between

(5 and (Er),,,

c%_ d

1,

150

time step Lo

Figure8 respect

Variation of percentagesaving to time step

in storage

requirement

with

2 CPl (mir

As mentioned previously, the truncation of the time integrals introduces an error. Let (Er),,,axbe the true error measure in terms of the maximum relative difference anywhere between the truncated and untruncated results of the problem; i.e.,

max

untruncated

We want

value -truncated

untruncated

to determine

value (13)

value

the relationship

between

6 and

Figure

Appl.

70

Math.

Determination

Modelling,

of rjopt

1987,

Vol. 11, December

415

Two-dimensional

transient

wave propagation

problems:

Summary and conclusions An efficient BEM for two-dimensional, transient, wave propagation problems has been outlined. The new approach is mainly aimed at improving the computation efficiency of a certain time-stepping scheme employed in BEM applications. We introduce a method for reducing computation time and storage requirements. The advantages of this approach are illustrated by two test cases where transient oscillation of a rectangular harbor is investigated. From the results of the numerical examples, we conclude that CPU time and storage requirements can be considerably reduced by truncating the time integrations. The negative effect of the truncation is noticeable only at the later time steps of the computation, where wave amplitudes tend to increase, although no phase shift occurs. A mathematical relationship for the determination of a truncation limit is derived in terms of a truncation guide parameter 6, defined through a typical function of the encountered problem. The correlation between 6 and the maximum error measure of the computed results, (Er),,,, indicates that an approximately linear relationship holds. Within the range 0.10 < 6 < 0.35, we found that (Er),,, < 6/2 for the problem studied. An optimal 6 exists through a balance of the requirement between error allowance and CPU time. In practice, a safe choice of this parameter may be taken to be 6 = 0.20, which should

416

Appt. Math. Modelling,

1987, Vol. 11, December

V. Demirel

and S. Wang

approximately correspond to a 7-8 % maximum error for most problems. If the leading wave of a transient wave system is important, the value of 6 may be extended up to 35%. As implied by Figure 7, a truncation of this magnitude may reduce the computation time more than 60% without affecting the accuracy of the leading wave. However, the trailing waves may have discrepancies of up to 15 %. It is anticipated that this simple but effective idea may be extended to other applications where a similar time-stepping scheme is employed. References Mansur, W. J. and Brebbia, C. A. Formulation of the boundary element method for transient problems governed by the scalar wave equation, Appl. Math. Modeiling, 1982, 6, 307 Mansur, W. J. and Brebbia, C. A. Numerical implementation of the boundary element method for two dimensional transient scalar wave propagation problems, Appl. Math. Modelling, 1982, 6, 299 Mansur, W. J. and Brebbia, C. A. Transient Elastodynamics, Topics in boundary element research, 1985, Chap. 5 Demirel, V. Transient linear long wave induced oscillations in arbitrary-shaped harbors: a boundary element approach, MS. Thesis, University of Miami, 1986 Pina, H. L. G. Time dependent potential problems, in Boundary element techniques in computer-aided engineering (ed C. A. Brebbia), Nato ASI. Series E, No. 84, Martinus Nijhoff, 1984 Demirel, V. and Wang, S. A numerical solution technique to transient linear long wave induced harbor oscillations using boundary element technique, Proc. 20th Int. Conf. on Coastal Engineering, Taipei, Taiwan, November, 1986