A new representation for the free-surface channel Green’s function

APOR 361

Applied Ocean Research 21 (1999) 17–25

A new representation for the free-surface channel Green’s function C.M. Linton Department of Mathematical Sciences, Loughborough University, Loughborough, Leicestershire LE11 3TU, UK Received 6 June 1998; accepted 18 September 1998

Abstract A new representation is derived for the free-surface Green’s function, appropriate to water wave problems in channels of constant depth and width, which is rapidly convergent throughout the fluid domain. The new representation depends on two arbitrary positive parameters a ~ and by varying these parameters the convergence characteristics of this new representation can be altered. Letting a and a~ tend to zero and a; results in the known eigenfunction expansion. The results of computations are presented showing the accuracy and efficiency of the new representation. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: Water waves; Green’s functions; Channels

1. Introduction As is well known, many problems concerning the interaction of water waves and structures can be formulated in terms of an integral equation involving an appropriate Green’s function, G…x; x0 †, see for example [1, Section 7.4.2]. The Green’s function can be thought of as the potential at a field point x resulting from an oscillating source at a point x0. The numerical solution of such an integral equation will be straightforward provided methods for the efficient evaluation of G are available and in a series of papers, Newman [2–4], has developed useful algorithms for the computation of Green’s functions satisfying the free-surface boundary condition which arises in linear water wave theory. Newman’s approach is to divide the fluid domain into various subregions and to use a different representation for the Green’s function in each of these subregions. Thus for the Green’s function appropriate to three-dimensional wave propagation in an unbounded ocean of constant depth he utilizes an eigenfunction expansion in an outer region where the horizontal distance between the source and field points is sufficiently large, but in the complementary inner region an approximation in terms of economized Chebyshev polynomials is used. Recently however, Linton [5] has derived expressions for the free surface Green’s function (in both two and three dimensions) which converge rapidly throughout the fluid domain (excluding, of course, at the singularity). The idea behind the derivation of these new representations for G is that the Green’s function satisfies Poisson’s equation in which the right-hand side is a product of Dirac delta func-

tions, and solutions to Poisson’s equation can be formally written as integrals of solutions to appropriate heat conduction problems. These latter solutions can be written in two ways, one of which converges rapidly for small values of the time, t, and one which converges rapidly for large t. When it comes to integrating this solution, the range of integration, …0; ∞†; is split and the different expressions are used in the appropriate parts of the integral. The result is a one-parameter family of representations for G (the positive parameter a governing the point at which the integration is split) and by varying this parameter the convergence characteristics of the representation can be altered. This family of representations is related to the eigenfunction expansion in that the latter is recovered by letting a ! 0. By taking non-zero values of a however, rapidly convergent expressions can be obtained in situations where the eigenfunction expansion converges slowly, if at all. In this article we consider the case of a channel of constant depth and width. This is a situation of considerable practical importance since experiments to determine the hydrodynamic characteristics of structures are often performed in wave tanks though a knowledge of the behaviour in the open ocean is required. Thus it is important to determine quantitatively the effect of the channel walls. For some special geometries involving vertical cylinders mathematical progress can be made using multipole expansions [6] but for more realistic structures the integral equation approach based on the Green’s function can be used, as in Kashiwagi [7] who considered an infinitely deep channel. One way to obtain the Green’s function for such a situation is to take the Green’s function for a horizontally

0141-1187/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0141-118 7(98)00032-7

18

C.M. Linton / Applied Ocean Research 21 (1999) 17–25

unbounded fluid and to write a sum representing the images of this singularity in the channel walls. However the resulting series converges extremely slowly and is useless from a computational point of view. (This is not to say that the method of images has no role to play in the numerical solution of water wave channel problems, it was used to good effect in [8] for example). An eigenfunction expansion can also be obtained which is computationally efficient for some parameter values but which converges slowly if the distance along the channel between the source and field points is too small. An alternative representation, suitable in a different part of the fluid region, was derived in [9] but even taken together these representations do not provide attractive methods for the computation of G at all points in the fluid. A representation which does provide an efficient method for the computation of G throughout the fluid region can however be obtained using the method of [5] and in Section 2 we will show how this representation can be obtained directly from the eigenfunction expansion. As a part of this process the Green’s function for the two-dimensional Helmholtz equation in the region between two parallel walls is required and we will use a rapidly convergent representation for this function given in [10]. Results of computations showing the accuracy and efficiency of the method are given in Section 3. All the computations were performed using Mathematica.

72 GN ˆ d…x†d…y ⫺ h†d…z ⫺ z† 0 ⬍ h ⬍ b; ⫺h ⬍ z ⱕ 0;

2GN ˆ0 2z

…2:6† where V ⬅ f ⫺∞ ⬍ x ⬍ ∞; ⫺∞ ⬍ y ⬍ ∞; ⫺h ⬍ z ⱕ 0g; together with the free-surface and bottom boundary conditions, Eqs. (2.3), (2.4) and the same radiation condition as for GN. It is readily demonstrated that, with r ˆ ‰x2 ⫹ y2 ⫹ …z ⫺ z†2 Š1=2 ; G⬃⫺

1 4pr

on y ˆ 0; b

in VN ;

…2:1†

…2:2†

on z ˆ 0;

…2:3†

on z ˆ ⫺h;

…2:4†

as r ! 0;

…2:7†

provided z ⬍ 0, but that if z ˆ 0, G⬃⫺

1 2pr

as r ! 0:

…2:8†

This Green’s function can be written as an eigenfunction expansion, see [9] for example. Gˆ⫺

∞ X

Zm

∞ X e⫺lmn 兩x兩 1 cosan y 2lmn d n nˆ0

Zm

X e⫺lmn 兩x兩 Y; 2lmn n n

mˆ0 ∞ X mˆ0

P Throughout this article n will be used to represent the sum over all integers n and we will make extensive use of the Dirac delta function d (·). We consider the three-dimensional fluid domain VN ⬅ { ⫺ ∞ ⬍ x ⬍ ∞; 0 ⬍ y ⬍ b; ⫺h ⬍ z ⱕ 0}; with the undisturbed free surface being z ˆ 0 so that the Green’s function representing an oscillating point source at x ˆ 0; y ˆ h; z ˆ z is Re{GN e⫺ivt }; where GN is the solution to

2GN ˆ KG 2z

n

ˆ⫺

2. Derivation

2GN ˆ0 2y

where G…x; y; z; z† represents a source periodic in y, the period being 2b. Below we will work entirely with this periodic Green’s function as it is simpler than GN above, though for convenience we will label the periodicity as d. It then satisfies X 72 G ˆ d…x†d…z ⫺ z† d…y ⫺ nd† in V; ⫺h ⬍ z ⱕ 0;

…2:9†

where 10 ˆ 1; 1n ˆ 2; n ⱖ 1; Yn ˆ eian y =d;

an ˆ 2np=d;

Zm ˆ Nm⫺1 cos mm …z ⫹ h†cos mm …z ⫹ h†;   h sin 2mm h 1⫹ ; Nm ˆ 2 2mm h

…2:10†

…2:11†

m0 ˆ ⫺im where m is the positive root of m tanhmh ˆ K;

…2:12†

mm ; m ˆ 1; 2; …; are positive and satisfy the equation mm tan mm h ˆ ⫺K;

…2:13†

and

lmn ˆ …a2n ⫹ m2m †1=2 ; l0n ˆ …a2n ⫺ m2 †1=2 ˆ ⫺i…m2 ⫺ a2n †1=2 :

…2:14†

and we require GN to behave like outgoing waves as 兩x兩 ! ∞. Here K ˆ v 2/g where v is the angular frequency of the waves and g is the acceleration owing to gravity. This Green’s function can be written as

Only a finite number of terms contribute to the imaginary part of G, namely the pairs m, n for which m ˆ 0 and a2n ⬍ m2 : Thus if s is a nonnegative integer such that 2sp ⬍ md ⬍ 2…s ⫹ 1†p; we have s   Z X 1n 2 2 1=2 cos … m ⫺ a † 兩x兩 cos an y: Im G ˆ ⫺ 0 n 2d nˆ0 …m2 ⫺ a2n †1=2

GN …x; y; z; h; z† ˆ G…x; y ⫺ h; z; z† ⫹ G…x; y ⫹ h; z; z†; …2:5†

…2:15†

C.M. Linton / Applied Ocean Research 21 (1999) 17–25

Provided 兩x兩 is not too small, the series in (2.9) converge rapidly, but when x is close to zero it becomes computationally inefficient. When x ˆ 0, the double summation in (2.9) converges extremely slowly, but the Poisson summation formula, which can be formally written as X

einu ˆ 2p

n

X

d…u ⫹ 2np†;

…2:16†

n

can be used to produce a more computationally useful representation. Thus, for any positive c, 1 X eian y 1 X Z∞ eiu…y⫺nd† ˆ du 2d n …a2n ⫹ c2 †1=2 4p n ⫺ ∞ …u2 ⫹ c2 †1=2 ˆ

1 X K …c兩y ⫺ nd兩†; 2p n 0

…2:17†

where we have used an integral representation for the modified Bessel function K0(·) given in Ref. [11, Eq. 8.432(5)]. It follows that when x ˆ 0   Z0 X ian y 1 1 e ⫺ 2 Gˆ⫺ 2d n l0n …an ⫹ m2 †1=2 ⫺

∞ X Z0 X eian y Zm X e i an y ⫺ 2d n …a2n ⫹ m2 †1=2 2d n …a2n ⫹ m2m †1=2 mˆ1

…2:18†

ˆ⫺

  ∞ Z0 X 1 1 1n cos an y ⫺ 2 2d nˆ0 l0n …an ⫹ m2 †1=2

…2:19† The sums of Bessel functions converge rapidly since K0 …x† ⬃ …p=2x†1=2 exp…⫺x† as x ! ∞; and the terms in the first summation in (2.19) are O…n⫺3 †: An alternative representation for the free-surface channel Green’s function was derived in [9] where it was shown that, for z ⬍ 0, 2pdG ˆ ln

r ⫺ r1

0

f0 …t; z; z†cos xt dt ⫺ 2

∞ X

⫹

e⫺ht sinh zt sinh zt t

…2:22†

and

ÿ
fn …t† ˆ K ⫹ an t ean zt cosh an …z ⫹ h†t ÿ
⫹ e⫺an …z⫹h†t K sinh an zt ⫹ an t cosh an zt :

…2:23†

In (2.20) the contours of integration are indented below any poles on the real axis, but only a finite number of the integrals have such poles and so the imaginary part can again be written as a finite sum. Whereas the eigenfunction expansion (2.9) converges rapidly provided 兩x兩 is large enough, the representation (2.20) can be used provided r /d is not too small and the source and field point are not both close to the free surface. In contrast, the new representation for G derived below is rapidly convergent throughout the fluid domain. We will derive the new representation directly from the eigenfunction expansion for which the m ˆ 0 term is ÿ
⫺c x; y Z0; where

cˆ

X e⫺l0n 兩x兩 Y: 2l0n n n

…2:24†

This is simply the eigenfunctioin expansion for the Green’s function for the two-dimensional Helmholtz equation in periodic domains, i.e. c is the solution of   X ÿ
72xy ⫹ m2 c ˆ ⫺d… x† d y ⫺ nd : …2:25† This function can be represented in many different ways and is discussed in detail in [10] where it is shown that the following representation is the most useful from a computational point of view. With a~ being an arbitrary positive parameter we can write X X cˆ L0n Yn ⫹ Y~ n ; …2:26† n

L0n

cos…an y†

 fn …t; z; z†cos…an xg…t†† dt 1 g…t†…an t sinh an ht ⫺ K cosh an ht …2:20†  1=2 where g…t† ˆ t2 ⫺ 1 ; h i1=2 r1 ˆ x2 ⫹ …z ⫹ z†2 ; …2:21†

n

where

nˆ1

 Z∞ × K 0 … an r † ⫹

h i1=2 r ˆ x2 ⫹ …z ⫺ z†2 ;

2 cosh ht

19

ÿ
ÿ
cosh z ⫹ h t cosh z ⫹ h t t sinh ht ⫺ Kcosh ht #

n

∞ X Z X Zm X K0 …m兩y ⫺ nd兩† ⫺ K …m 兩y ⫺ nd兩†: ⫺ 0 2p n 2p n 0 m mˆ1

Z∞

f0 …t† ˆ

"

" ! ~ 1 l0n ad x l0n x ⫹ ˆ e erfc ~ 2 4l0n ad !# ~ l0n ad x ⫺l0n x ⫺ erfc ⫹e ~ 2 ad

…2:27†

and Za2 d2 =4 e⫺R2n =4t⫹m2 t dt 4pt 0 !   ∞ ~ 2p 1 X 1 mad R2n Ep⫹1 2 2 : ˆ 2 4p pˆ0 p! a~ d

Y~ n ˆ

…2:28†

20

C.M. Linton / Applied Ocean Research 21 (1999) 17–25

Here Rn ˆ ‰x2 ⫹ …y ⫺ nd†2 Š1=2 ;

and so

2 Z∞ ⫺t2 e⫺z e dt ⬃ p1=2 z p1=2 z 2

erfc…z† ˆ

G ˆ ⫺ cZ 0 ⫺

mˆ1

as z ! ∞; 兩arg z兩 ⬍ 3p=4

⫺

is the complementary error function and Z∞ e⫺x as x ! ∞ u⫺n e⫺xu du ⬃ E n … x† ˆ x 1

∞ X X

Yn Zm

n

mˆ1

Z∞ e⫺x2 =4t⫺l2mn t dt: …4pt†1=2 0

Lmn ˆ

ˆ

Z∞ a2 h2 =4

=4t⫺l2mn t

e⫺x dt …4pt†1=2

1 2p1=2 lmn

…2:30†

where the series in terms of the incomplete Gamma function is derived by expanding exp( ⫺ x 2/4t) in a power series and then ! ∞ X Za2 h2 =4 e⫺x2 =4t⫺l2mn t X Yn Zm Lmn ⫹ dt G ˆ ⫺cZ0 ⫺ …4pt†1=2 0 mˆ1 n …2:32†

mˆ1

⫺

∞ X mˆ1

n

Zm

Za2 h2 =4 e⫺x Xe 1=2 d …4pt† 0 n =4t⫺m2m t

wˆ

∞ X

e⫺mm t Zm : 2

…2:36†

mˆ0

This series converges rapidly for large values of t, but the integral in (2.35) is over values of t near the origin so we require a representation of w which converges rapidly for small t. Such an expansion was derived in [5], in which it was shown that ∞ 4 X X ÿ
e⫺…z⫺z† =4t e⫺…2h⫹z⫹z† =4t ⫹ ⫹ …⫺1†m Im xm;i : wˆ 1=2 1=2 …4pt† …4pt† mˆ1 iˆ1 2

2

…2:37†

xm;1 ˆ 2…m ⫺ 1†h ⫺ z ⫺ z;

…2:38†

xm;2 ˆ 2mh ⫺ z ⫹ z;

…2:39†

xm;3 ˆ 2mh ⫹ z ⫺ z;

…2:40†

xm;4 ˆ 2…m ⫹ 1†h ⫹ z ⫹ z;

…2:41†

and it can be seen that since –h ⬍ z ⱕ 0 and –h ⬍ z ⱕ 0, all the x m,i that occur in (2.37) are positive. In fact, for m ⱖ 1, 2…m ⫺ 1†h ⱕ xm;i ⬍ 2…m ⫹ 1†h:

ian y⫺a2n t

dt:

…2:33†

m ÿ
X ÿ
m! ÿ
…⫺1†j …2K †m⫺j I~m⫺j x Im x ˆ j! m ⫺ j ! jˆ0

…2:43†

…2:44†

with the initial values ÿ
e⫺x =4t ; I~0 x ˆ …4pt†1=2

…2:45†

  ÿ
1 2 x 1=2 ⫺ Kt : I~1 x ˆ ⫺ eK t⫺Kx erfc 2 2t1=2

…2:46†

2

Now from the Poisson summation formula we can derive the identity X ⫺ y⫺nd 2 =4t 1 X ⫺a2n t ian y 1 † ; e e ˆ e … 1=2 d n …4pt† n

…2:42†

The functions Im(x ) that occur in (2.37) are defined by the formula

and the recurrence relation
ÿ
ÿ
ÿ ÿ
…n ⫺ 1†I~n x ˆ 2tI~n⫺2 x ⫹ x ⫺ 2Kt I~n⫺1 x ;

Lmn Yn Zm 2

…2:35†

The quantities x m,i are defined by

!   ∞ X …⫺1†p xlmn 2p 1 l2mn a2 h2 ⫺ p; ; G 2 p! 2 4 pˆ0

∞ X X

Za2 h2 =4 e⫺R2n =4t⫺m2m t dt: Zm 4pt 0

…2:29†

…2:31†

ˆ ⫺cZ0 ⫺

n

Lmn Yn Zm

n

Next we consider the function

The key to the derivation of the new representation for G is to split the integral in this expression at some arbitrary point depending on a positive parameter a and treat the two resulting integrals differently. Thus for m ⱖ 1 we write 2

∞ X X mˆ1

is the exponential integral. The series in (2.28) is readily obtained from the integral by expanding exp(m 2t) in a power series. The parameter a~ controls the respective rates of convergence of the two summations in (2.26): increasing a~ makes the first summation converge faster and the second slower. If we let a~ ! 0 in (2.26) we recover the eigenfunction expansion (2.24). If we separate the m ˆ 0 term from the rest of the eigenfunction expansion for G and use [11, Eq. 3.325], we can write G ˆ ⫺cZ0 ⫺

∞ X X

…2:34†

C.M. Linton / Applied Ocean Research 21 (1999) 17–25 Table 1 Bound on 兩L2n 兩 from (3.2), Kh ˆ 1

and

a

Lmn ˆ 3.9 × 5.1 × 2.6 × 2.0 × 3.7 × 9.9 × 3.4 × 3.9 × 2.3 × 9.2 ×

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

⫺177

10 10 ⫺46 10 ⫺21 10 ⫺12 10 ⫺8 10 ⫺6 10 ⫺4 10 ⫺3 10 ⫺2 10 ⫺2

4 Za2 h2 =4 e⫺R2n =4t X ÿ
Im xm;i dt: 4pt iˆ1 0

21

…2:52†

The imaginary part of G comes entirely from L0n, given by (2.27), and it is only the finite number of values of n for which l 0n is imaginary which contribute, again leading to (2.15). If we let a and a~ tend to zero in (2.50) we recover the eigenfunction expansion (2.9). Note that

Yn ⫺ Y~ n ˆ

Za2 h2 =4 e⫺R2n =4t⫹m2 t dt 4pt a~ 2 d 2 =4

…2:53†

which vanishes if we choose a˜ ˆ ah/d.

In particular ÿ
e⫺x =4t ; I0 x ˆ …4pt†1=2 2

3. Results

2   ÿ
e⫺x =4t x K 2 t⫺K x 1=2 ⫺ Ke erfc ⫺ Kt : I1 x ˆ ⫺ 2t1=2 …4pt†1=2

…2:47†

In terms of w we have G ˆ ⫺cZ0 ⫺

∞ X X mˆ1

⫺

Lmn Yn Zm

n

2  X Za2 h2 =4 e⫺Rn =4t  2 w ⫺ em t Z0 dt; 4pt 0 n

…2:48†

and if we define Za2 h2 =4 e⫺R2n =4t⫹m2 t dt Yn ˆ 4pt 0 !   ∞ 1 X 1 mah 2p R2n Ep⫹1 2 2 ; ˆ 4p pˆ0 p! 2 a h

…2:49†

and substitute for c from (2.26) and for w from (2.37), we obtain our new representation for the free-surface channel Green’s function: Gˆ

When it comes to the numerical evaluation of (2.50), a major simplification arises from the fact that, provided we take a sufficiently small, all the terms Lmn, m ⱖ 2 can be neglected. These terms decay rapidly as m increases and it is not difficult to show that provided Kh ⬍ 2a ⫺2,

∞ X X Xÿ
Yn ⫺ Y~ n Z0 ⫺ Lmn Yn Zm n

⫺

mˆ0

X

"

n

⫹

∞ X

1 r erfc n ah 4prn # !m ⫺1

!

n

1 r 0n ⫹ erfc ah 4pr 0n

!

兩L2n 兩 ⬍

Lmn ;

…2:50†

mˆ1

where h i1=2 ÿ
rn ˆ x2 ⫹ y ⫺ nd 2 ⫹…z ⫺ z†2 ; r 0n

h

ÿ


2

i 2 1=2

ˆ x ⫹ y ⫺ nd ⫹…2h ⫹ z ⫹ z† 2

In this section we will demonstrate the accuracy and efficiency of the formula (2.50) for the free-surface channel Green’s function, G. For computational purposes we need to determine the optimal values for the parameters a and a˜, and truncate the various infinite series that occur after an appropriate number of terms. The relative merits of using a particular method for evaluating G in a particular instance will of course depend on the precise numerical implementation and the availability or otherwise of routines to evaluate any necessary special functions. Here we will only examine the number of terms required to achieve a desired accuracy; the actual CPU time taken to evaluate the individual terms will not be addressed except by a few general remarks. In the calculations we evaluated the integrals Yn ⫺ Y~ n ; and Lmn ; m ⱖ 1, using standard numerical quadrature routines, rather than from the alternative series representations, which is straightforward because of the nature of the integrands. There is one exception to this and that is the evaluation of Lmn when x ˆ 0, in which case   1 l ah Lmn ˆ erfc mn : …3:1† 2lmn 2

…2:51† ;

   a2 ⫺4=a2  e 2 ⫹ a2 8Kh ⫹ 49K 2 h2 ⫹ 4K 3 h3 a4 ; 4p …3:2†

for any values of x, y, z,z and n. An example of the numerical value of the right-hand side of this inequality as a function of a is shown in Table 1, in which we have taken Kh ˆ 1. Note that this bound is fairly crude and that if any of x, z, z are not close to zero, or y is not close to zero or d, or n ⬎ 1, the bound will reduce significantly. If

22

C.M. Linton / Applied Ocean Research 21 (1999) 17–25

Table 2 Truncation parameters required to achieve on accuracy of 8 decimal places in Re G using (3.3), when y/h ˆ 0.5, z/h ˆ ⫺ 0.5, z=h ˆ ⫺0:2, Kh ˆ 1 and d/h ˆ 1 M, N1, N2 a

x/h ˆ 0.1

x/h ˆ 0.01

x/h ˆ 0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

40,9,— 12,6,0 10,5,1 6,3,1 5,2,1 4,2,2 4,2,2 3,1,3

350,25,— 14,8,0 8,5,1 6,3,1 5,2,2 4,2,2 4,2,2 2,1,3

10, 110, 15 (from Eq. (3.7)) 10,10,0 8,5,1 6,3,1 5,2,1 4,2,2 4,2,2 2,1,3

Kh is increased however, the bound will also increase, though not particularly rapidly. Numerical calculations show that for values of Kh in the range 0 ⱕ Kh ⱕ 10 (which is likely to include any value of practical interest) a choice of a ⱕ 0:4 will ensure that the first 8 decimal places (at least) of G will be unaffected by neglecting the terms Lmn, m ⱖ 2. We now introduce the truncation parameters M, N1 and N2 so that, neglecting the aforementioned terms, we can write N1 M X 1 X Zm 1 L cos an y d mˆ0 nˆ0 n mn

Gˆ⫺

N2 X

⫹

(



4  X Yn ⫺ Y~ n Z0 ⫹ Qn ⫺

nˆ⫺ N2

iˆ⫺ 1

 ) Rn;i 1 erfc ah 4pRn;i

⫹… …3:3† where  1=2 Rn;⫺1 ˆ R2n ⫹ …z ⫺ z†2 ; 

Rn;0 ˆ R2n ⫹ …2h ⫹ z ⫹ z†2  1=2 Rn;i ˆ R2n ⫹ x21;i ;

1=2

…3:4† ;

i ˆ 1; …; 4

…3:5†

and Qn ˆ ⫺

×

K Zah=2 K 2 u2 e ⫺ R2n =4u2 2p 0 4 X iˆ1

e⫺K x1;i erfc



 x1;i du ⫺ Ku : 2u u

…3:6†

When a ˆ 0 Eq. (3.3) reduces to a truncated version of the eigenfunction expansion, (2.9), and N2 ceases to be relevant. When x ˆ 0 we also have the representation ! N1 Z0 X 1 1 1 cos an y ⫺ ÿ Gˆ⫺
2d nˆ0 n l0n a2n ⫹ m2 1=2 1 ⫺ 2p

N2 X

(

) M
X ÿ
Z0 K0 m兩y ⫺ nd兩 ⫹ Zm K0 mm 兩y ⫺ nd兩 : ÿ

nˆ⫺ N2

mˆ1

…3:7† In the computations below we will determine the truncation parameters required (where the truncation parameters are large they have been rounded up to the nearest 10) to achieve an accuracy of 8 decimal places in the real part of G (the imaginary part is determined by a finite number of terms). This is a greater accuracy than is usually required in engineering practice but validation of programs and numerical methods requires more accurate calculations and also, as the derivative of G is often required, the truncation parameters calculated below will be suitable in any practical situation for the computation of the first spatial derivatives of G using term by term differentiation. (The effect of differentiation is to introduce algebraically growing factors into the exponentially decaying series; typically one can expect one or two fewer decimal places of accuracy from the same truncation parameters.) ~ which Formula (3.3) is valid for all positive values of a, appears through the terms L0n and Y~ n . If we analyse the behaviour of these terms as 兩n兩 ! ∞ we find that    L0n ˆ O n⫺2 exp ⫺n2 p2 a~2 ⫹ 2np兩x兩=d ; …3:8†    Y~ n ˆ O n⫺2 exp ⫺n2 =a~2 :

…3:9†

Thus increasing a~ will make Sn L0n converge faster and Sn Y~ n slower; when x ˆ 0 the rates of convergence balance at a~ ˆ p⫺1=2 ⬇ 0:56. This suggests that a value of a~ in the vicinity of 0.5 is desirable. Computations from [5] suggest that a similar value is desirable for a, though in this case the neglect of the terms Lmn, m ⱖ 2 means that we must keep a less than 0.5. Finally we note that we would like to choose a~ ˆ ah=d so that Yn ⫺ Y~ n vanishes, and this is clearly sensible if d/h ⬇ 1, but for large or small aspect ratios the need to compute Yn ⫺ Y~ n must be balanced against using non-opti~ Note that the width of the channel in mal values for a or a. (2.1)–(2.4) is actually b ˆ d/2, hence d/h ˆ 1 corresponds to a channel of width equal to half the depth. We will begin by considering the most favourable aspect ratio, d/h ˆ 1, and will set a~ ˆ a. We will choose y/h ˆ 0.5, z/h ˆ ⫺ 0.5, z=h ˆ ⫺0:2 and Kh ˆ 1, values which cause no particular problems, and investigate the convergence properties of the representations described before, for various values of a and x/h. The results are displayed in Table 2. For the three values of x/h shown, the same accuracy can be achieved by truncating the representation (2.20) at n ˆ 8.

C.M. Linton / Applied Ocean Research 21 (1999) 17–25 Table 3 Truncation parameters required to achieve on accuracy of 8 decimal places in Re G using (3.3), when x/h ˆ 0, y/h ˆ 0.5, z/h ˆ z /h ˆ and d/h ˆ 1 M, N1, N2 a

Kh ˆ 0.5

0 0.1 0.2 0.3 0.4 0.5 a

12,80,22 11,12,0 10,6,1 7,4,1 6,3,2 4,2,2

a

Kh ˆ 1 12,200,14 10,12,0 9,6,1 7,4,1 5,3,2 4,2,2

Kh ˆ 4 a

9,300,8 a 10,11,0 9,6,1 7,4,1 5,3,2 4,2,2

From Eq. (3.7).

The eigenfunction expansion converges rapidly for sufficiently large 兩x兩 and so when 兩x兩/h is large it is sensible to set a ˆ 0 in Eq. (3.3). The precise meaning of ‘large’ depends on the numerical implementation but in the case considered here it was found that the computation of G was quicker using the eigenfunction expansion when x/h ⱖ 0.1 but quicker using a non-zero value of a when x/h ⱕ 0.01, the transition occurring somewhere in-between. As x/h decreases towards zero the eigenfunction expansion converges too slowly and when x ˆ 0 thousands of terms are required to obtain 8 decimal place accuracy. However, when x ˆ 0 we have the alternative, Eq. (3.7), and the truncation parameters for this representation are given in Table 2. The results in Table 2 clearly demonstrate how increasing a reduces the number of terms required in the eigenfunction part of the representation for G (M and N1) at the expense of the need to calculate terms in the second summation (N2). What is remarkable, given the change in the convergence of the eigenfunction expansion as x/h ! 0, is the lack of sensitivity of the required truncation parameters to the value of x/h. As will be seen throughout the results that follow, varying the physical parameters has very little effect on the convergence characteristics of the new representation and this is an attractive feature. It was suggested above that optimal performance would be achieved near a ˆ 0.5. For the parameters used in Table 2

Table 4 Truncation parameters required to achieve an accuracy of 8 decimal places in Re G using (3.3), when x/h ˆ 0, y/h ˆ 0.1, z/h ˆ ⫺ 0.5, z /h ˆ 0.2, Kh ˆ 1 and d/h ˆ 0.2 M, N1, N2 a~

a ˆ 0.1

a ˆ 0.2

a ˆ 0.3

a ˆ 0.4

a ˆ 0.5

0.5 1.0 1.5 2.0 2.5

22,1,1 22,2,4 22,2,6 22,2,7 22,2,9

12,2,4 12,1,4 12,1,6 12,1,8 12,1,10

8,2,5 8,1,5 8,0,6 8,0,8 8,0,10

6,2,7 6,1,7 6,0,7 6,0,8 6,0,10

4,2,9 4,1,9 4,0,9 4,0,9 4,0,10

23

the restriction a ⬍ 0.5 is not necessary and in the table we have included data for a ˆ 0.6 and a ˆ 0.7 (8 decimal place accuracy was not possible with a ˆ 0.8). It is clear that 0.5 is an underestimate for the optimal value, though any a between 0.3 and 0.7 gives a fast and efficient method for the accurate evaluation of G. For the remainder of the results we will concentrate on x ˆ 0 since it is when 兩x兩/h is small that the new representation is most useful. If we use x ˆ 0 rather than a small but non-zero value of 兩x兩/h we can compare our new representation, (3.3), with the form (3.7). Moreover, Table 2 shows that the convergence characteristics of (3.3) are unaffected by whether x is exactly zero or not, so the conclusions reached below for x ˆ 0 are equally appropriate for any small value of 兩x兩/h. As has already been mentioned, the neglect of the terms Lmn, m ⱖ 2, has its greatest effect when z ˆ z ˆ 0, i.e. when both the source and the field point lie in the free surface, and in such circumstances the representation (2.20) is not valid. The results in Table 3 show that the new representation still converges rapidly in this case. In this example we need a ⱕ 0.5 to ensure an accuracy to 8 decimal places when Kh ⱕ 4. The table shows results for three different values of the nondimensional frequency parameter covering the range of greatest practical interest and it is apparent that the convergence characteristics of the new representation are not sensitive of the value of Kh in this range. Truncation parameters for (2.19) are also shown and the time taken to compute G from (2.19) was roughly comparable with using the new representation with a ˆ 0.5. So for channels with d/h near unity, i.e. b ⬇ h/2, one should choose a~ ˆ ah=d and a ⬇ 0.5 for optimal performance. What about channels for which the ratio d/h is large or small? We will begin with the case of a narrow channel and take d/h ˆ 0.2. In this case, in order that Yn ⫺ Y~ n vanishes we need a~ ˆ 5a which may cause problems. However, the value of a is more important than the value of a~ since it controls the convergence of the more computationally expensive part of the representation (3.3) and is also restricted because of the neglect of the terms Lmn, m ⱖ 2. Table 4, in which y/h ˆ 0.1, z/h ˆ ⫺ 0.5, z /h ˆ ⫺ 0.2 and Kh ˆ 1, shows the truncation parameters required to achieve an accuracy of 8 decimal places for various values of a and ~ The value of M does not depend on a~ and we see that this a. is smallest when a ˆ 0.5. The case a~ ˆ 5a consists of the diagonal entries in the table and these are seen to be close to optimal. It is worth noting however that the actual time taken to evaluate Yn ⫺ Y~ n when a~ 苷 ah=d is very small, and choosing a~ ˆ ah=d may well be less important than simply trying to restrict the number of terms required. For the values of the parameters used in Table 4 the truncation parameters required in the representation (3.7) are (M, N1, N2) ˆ (60,30,120), which makes that method unattractive in this case, and the representation (2.20) converges to 8 decimal places with n ˆ 2. In fact (2.20) converges to 7 decimal places with n ˆ 1 and for narrow

24

C.M. Linton / Applied Ocean Research 21 (1999) 17–25

Table 5 Truncation parameters required to achieve an accuracy of 8 decimal places in Re G using (3.3), when x/h ˆ 0, y/h ˆ 2.5, z/h ˆ ⫺0.5, z /h ˆ ⫺0.2, Kh ˆ 1 and d/h ˆ 5 M, N1,N2 a~

a ˆ 0.1

a ˆ 0.2

a ˆ 0.3

a ˆ 0.4

a ˆ 0.5

0.05 0.1 0.2 0.5

2,51,— 2,51,— 2,51,1 2,51,2

2,26,— 2,26,— 2,26,1 2,26,2

2,22,— 2,17,— 2,17,1 2,17,2

2,22,— 2,13,— 2,13,1 2,13,2

2,22,— 2,11,— 2,10,1 2,10,2

channels (and 兩z兩/h and 兩z 兩/h not both small) this method becomes competitive. Results for a wide channel with d/h ˆ 5 are shown in Table 5 in which a dash for N2 indicates that no terms are required in the second summation. In this case the required value of M, which only depends on a, is much less than for the narrow channel above, but the number of cross-channel modes required (governed by N2) is much larger. Problems can arise if a~ is taken too small, but with a ˆ 0.5, a~ ˆ ah/d ˆ 0.1 is both the desired and the best choice. The representation (3.7) requires the truncation parameters (M,N1,N2) ˆ (2, 250,3) for these parameter values and the representation (2.20) converges too slowly to be of any use. Now G is singular at r ˆ 0 and the nature of the singularity is described by (2.8) or (2.7) depending on whether z is zero or not. In applications it is important to be able to accurately determine the regular part of G, since the singular part can be treated analytically. Hence it is usually desirable to have a representation which contains the singularity explicitly. None of the representations above have this property however, but Table 6 shows that the representation (3.3) is sufficiently accurate and computationally efficient for this not to matter. In the table we show the value of the regular part of Re G to 8 decimal places and the truncation parameters required to achieve this accuracy, as the field point approaches the singular point along the y-axis. Two sets of parameter values are considered: in both cases x ˆ 0; Kh ˆ d=h ˆ 1 and a~ ˆ a ˆ 0:4 and the values of z and z are chosen to make the accurate evaluation of G as awkward as possible. In the example on the left z=h ˆ z=h ˆ ⫺0:01 so that G ⬃ ⫺1=4py as y ! 0 whereas in the example on the right z=h ˆ z=h ˆ 0 so that G ⬃ ⫺1=2py as y ! 0. The appropriate behaviour is clearly demonstrated, as is the fact that the number of terms required to achieve an

accuracy of 8 decimal places in the regular part of G does not increase as we approach the singular point. Finally in this section we will list some simple rules by which sensible (though not necessarily optimal) values of the parameters a and a~ can be chosen. First, if 兩x兩/h ⱖ 0.05 use a ˆ a~ ˆ 0, which is just the eigenfunction expansion. For 0 ⱕ 兩x兩=h ⬍ 0:05 and Kh ⱕ 4 and channels with ‘moderate’ aspect ratios (i.e. channels for which the width b and depth h satisfy 0.1 ⱕ b/h ⱕ 2.5), use a ˆ 0.5 and a~ ˆ ah=d. For ‘narrow’ channels use a~ ˆ 2 rather than a~ ˆ ah=d and for ‘wide’ channels use a~ ˆ 0:2. When Kh ⬎ 4 use a ˆ 0.4 rather than 0.5.

4. Conclusion We have derived a new representation for the free-surface Green’s function suitable for problems set in a channel of constant depth and width in the form of a one-parameter family of formulas depending on two positive parameters ~ The convergence characteristics of the new reprea and a. ~ are sentation (which can be altered by varying a and a) excellent and what is more, largely insensitive to changes in the underlying physical parameters. The well-known eigenfunction expansion is a special case of the new representation, which can be obtained by letting a and a~ tend to zero, and is the best representation from which to compute G when the distance down the channel between the source and field point is sufficiently large. The great advantage of the new representation is in the region where this distance is small so that the eigenfunction expansion converges too slowly to be useful. Moreover, when the source and field points both lie close to the free surface so that the alternative representation derived in [9] cannot be used, the new representation still converges rapidly. The new representation performs best when the channel width is about half the depth, but is still computationally efficient for large and small channel aspect ratios. When the source and field points are at the same station along the channel the new representation is particularly efficient and the eigenfunction expansion converges too slowly to be of any use. However for this case the eigenfunction expansion can be manipulated, using standard techniques to produce an accelerated series which is particularly useful for narrow channels.

Table 6 Truncation parameters required to achieve an accuracy of 8 decimal places in the regular part of G using (3.3), when x/h ˆ 0, Kh ˆ d/h ˆ 1 and a ˆ a~ ˆ 0.4

y/h

z=h ˆ z=h ˆ ⫺0:01 M,N1,N2

0.1 0.01 0.001 0.0001

5,2,1 6,3,1 5,3,1 6,3,1

Re G ⫹ (4p y) ⫺1 ⫺ 0.49316979 ⫺ 3.44206061 ⫺ 3.86630525 ⫺ 3.87131946

z=h ˆ z=h ˆ 0 M,N1,N2

Re G ⫹ (2p y) ⫺1

5,3,1 5,3,1 5,3,1 5,3,1

0.26096468 ⫺ 0.11267629 ⫺ 0.48051461 ⫺ 0.84712512

C.M. Linton / Applied Ocean Research 21 (1999) 17–25

References [1] Mei CC. The applied dynamics of ocean surface waves. New York: Wiley-Interscience, 1983. [2] Newman JN. Algorithms for the free-surface Green function. J Engng Math 1985;19:57–67. [3] Newman JN. The evaluation of free-surface Green functions. In 4th Intl Conf. On Numerical Ship Hydrodynamics, Washington, DC, 1985. [4] Newman JN. Approximation of free-surface Green functions. In: Martin PA, Wickham GR, editors. Wave Asymptotics, Cambridge University Press, 1992. pp. 107–135. [5] Linton CM. Rapidly convergent representations for Green’s functions for Laplace’s equation. Proc. Roy. Soc. Lond. A (to appear). [6] Linton CM. The use of multipoles in channel problems. Chapter 2.

[7]

[8] [9] [10] [11]

25

In: Mandal BN, editor. Mathematical Techniques for Water Waves, International Series on Advances in Fluid Mechanics. Southampton: Computational Mechanics Publications, 1997, pp. 45–78. Kashiwagi M. Radiation and diffraction forces acting on an offshorestructure model in a towing tank. Intl Journal of Offshore and Polar Engng 1991;1:101–107. Yeung RW, Sphaier SH. Wave-interference effects on a truncated cylinder in a channel. J. Engng. Math. 1989;23:95–117. Linton CM. On the free-surface Green’s function for channel problems. Applied Ocean Research 1993;15:263–267. Linton CM. The Green’s function for the two-dimensional Helmholtz equation in periodic domains. J Engng Math 1998;33:377–402. Gradshteyn IS, Ryzhik IM. Tables of Integrals, Series and Products. New York: Academic Press, 1965.