Mathematical models for simple harmonic linear water waves: Wave diffraction and refraction

Mathematicar Models for

Simple Harmonic Linear Water Waves

^lZ..'y^)2j

This thesis is also published as Delft Hydraulics Laboratory

Publication No. 163 (1976).

IVIathematical IVlodels

for

Simple Harmonie Linear Water Waves

Wave Diffraction and Refraction

P R O E F S C H R I F T ter verkrijging van

de graad van doctor in de

technische wetenschappen

aan de Technische Hogeschool Delft

op gezag van de rector magnificus

Prof. dr. ir. H. van Bekkum,

hoogleraar in de afdeling der

Scheikundige Technologie,

voor een commissie aangevi/ezen

door het college van dekanen

te verdedigen op

woensdag 7 april 1976 te 16.00 uur door

J U R I C O R N E L I S W I L L E M B E R K H O F F

wiskundig-ingenieur

r

Dit proefschrift is goedgekeurd door de promotor

PROF. DR. E. VAN S P I E G E L

This thesis was prepared while the author was employed at the

De Voorst Laboratory of the Delft Hydraulics Laboratory. The

investigation reported herein was partly conducted under contract

for the Dutch Ministry of Public Works (Rijkswaterstaat)

Contents

I INTRODUCTION 1

1.1 Linear Surface Waves and Mathematical Models '

1 .2 Literature 3

1.3 Outline

<*

2

MATHEMATICAL FORMULATION 5

2.1 General Equations 5

2.2 Energy Flux 8

2.3 Potential Formulation 9

2.4 Boundary Conditions 10

2.5 Linear Equations 1'

2.6 Harmonic Solutions '2

2.7 Fundamental Solution 1^

2.8 Solution with Singularity Distributions 16

3 DIFFRACTION 20

3.1 The Diffraction Equation 20

3.2 Boundary Conditions 22

3.3 Some Analytical Solutions 25

3.4 Wave Penetration into Harbours 28

3.5 Harbour Oscillation 31

3.6 Numerical Method of Solution 34

3.7 Some Results 37

3.8 Discussion 38

4 REFRACTION 45

4.1 The Refraction Equations 45

4.2 Boundary Conditions 54

4.3 Some Analytical Solutions 56

4.4 Numerical Method of Solution 58

4.5 Some Results 60

4 . 6 Discussion 64

r

_{5 REFRACTION-DIFFRACTION 65} 5.1 The Refraction-Diffraction Equation 65

5.2 Boundary Conditions 69 5.3 Some Analytical Solutions 69 5.4 Variational Formulation 76 5.5 Numerical Method of Solution 79

5.6 Some Results 85 5.7 Discussion 96 CONCLUSIONS REFERENCES NOTATION SUMMARY SAMENVATTING

1 Introduction

1.1 Linear Surface Waves and Mathematical Models

Propagation of water waves is a phenomenon that interests oceanographic, meteoro-logie and civil engineers.

The phenomenon in all its aspects is a very complex one, because several physical processes play a part, such as the irregularity of the waves in nature, the break-ing of the waves and other non-linear effects, energy dissipation due to friction or turbulence, partial reflection, and the dependence of the wave celerity and direction of propagation on the bottom lay-out.

It is nearly impossible to give a mathematical description of the phenomenon of wave propagation in which all effects are included, and therefore some simplifica-tions must be made depending on the aspects which are of interest.

The applicability of a mathematical description often depends on the number of space dimensions involved with the problem. In the case of a one-dimensional model it is possible to include some non-linear effects without introducing too many difficulties into the method of solution. However, in the case of two-and three-dimensional models the mathematical formulation is often restricted to the linear simple harmonic wave theory.

The irregularity of the wave field is described then as a linear stochastic phenomenon, characterised by an energy density function.

For practical applications even a three-dimensional linear wave model is not applicable for problems dealing with wave propagation over a large area. A dimin-ishing of the number of space variables to the two horizontal coordinates would give a more applicable mathematical model. It is possible to formulate several two-dimensional models, depending on some assumptions about the form of the solu-tion domain and the behaviour of the solusolu-tions.

This present study deals with three mathematical models for wave propagation problems in the horizontal plane. These models are called diffraction, refraction and the combined form of refraction and diffraction. With the aid of Figure 1.1 it is possible to show the restrictions and assumptions that lead to the several mathematical models. It is an illustration of the characteristic quantities and parameters, which are important in wave propagation problems.

- ^ > ~ ^

S = H / I O M' D/IO a. = H/D £ = ff/M ffD

-->

D = mean water depth H = wave height •fo= 9 / ' ^ '

cr = mean slope of the bottom over a distance D

Figure 1.1 Characteristic Lengths and Parameters

For a linear wave model it is assumed that the wave steepness y ~ ^/l is small and in the case of a small water depth the ratio a = H/D must also be small.

For the diffraction model it is assumed that the bottom slope a over a distance that corresponds to the water depth is so small, or the water is so deep,that there is no influence of the bottom slope upon the surface waves. This assumption can be described with the requirement e « 1 in which £ = ol\x whith y = D/1 .

For greater values of the parameter E there is some influence of the bottom va-riation upon the surface waves. In the refraction theory it is assumed that the variation of the bottom is moderately small and noticeable for the surface waves ( E < 1), Also it is assumed that the variation of the wave amplitude in the horizontal plane depends only on the bottom variation and can be neglected in first approxi-mation. This last assumption about the behaviour of the solution is a rigid one, which gives rise to the problem of a caustic, where the solution of the refrac-tion equarefrac-tions is not well defined. When caustics are present it is clear that neglecting the curvature of the wave amplitude in the horizontal plane is not allowed. Caustics can occur also in those cases when the slope of the bottom is very small, so they are not only caused by the bottom variation but also by the boundary conditions of the refraction equations.

In the combined refraction-diffraction mathematical model it is assumed again that the variation of the bottom is moderately small and noticeable for the surface waves (e < 1). However, the curvature of the wave amplitude is now taken

into account in some way, A two-dimensional equation in the horizontal plane is obtained by integration over the water depth.

It is the purpose of this study to give a survey of the mathematical models based on the simple harmonic linear wave theory with a derivation of the basic differential equations and with a description of the methods of solution.

1 .2 Literature

A general mathematical description of the problems in fluid dynamics dealing with wave propagation was given a long time ago and can be found in several hand-books such as Lamb (1963), Wehausen (1960), Stoker (1957) and many others. There are many workers in the field of linear wave theory. A literature

survey about diffraction of water waves can be found in an article of Newman (1972).

For three-dimensional linear wave problems, which are of interest for the com-putation of wave forces on structures, reference should be made to the work of We-hausen (1960), John (1950), Garrison (1972) and van Oortmersen (1972). Two-dimen-sional diffraction problems, such as wave penetration into harbours with a con-stant water depth and harbour oscillations, have been studied by Biesel and Ran-son (1961), Daubert and Lebreton (1965), Barailler and Gaillard (1967), Lee (1969), Lee and Raichlen (1971) and the author (1969).

The describing differential equations of the refraction problem usually are de-rived by means of the geometric laws of Snellius and on the assumption that the wave energy flux between two wave rays is constant. Publications in this field have been given by Munk and Arthur (1952), Griswold (1963), Lepetit (1964), Wil-son (1966), DobWil-son (1967), Skovgaard (1973) and many others. However, it is pos-sible to derive these equations by means of an asymptotic approximation, concern-ing which reference can be made to articles by Keller (1958), Battjes (1968) and the author (1970).

In the past many attempts have been made to derive a two-dimensional equation which describes the combined refraction-diffraction effects in an appropriate way. Pierson (1951), Eckart (1952), Battjes (1968), Biesel (1972) and Ito (1972) have made some suggestions, but none of these equations is applicable in the whole field from deep to shallow water. The author (1972) and, independently of him,Schönfeld (1972) have succeeded in deriving such an equation.

1.3 Outline

It is important in the use of mathematical models to know the restrictions and simplifications which have been made to obtain these models. It is also conven-ient when all these simplifications are mentioned in one study and therefore Chapter 2 starts with a derivation of the general equations of motion for an ideal fluid by means of the divergence theorem of Gauss. The Chapter continues with the required simplifications to come to the three-dimensional linear simple harmonic wave equations and gives a general method of solution based on

singular-ity distributions.

The two-dimensional diffraction, refraction and refraction-diffraction models are described in the Chapters 3, 4 and 5 respectively. In each Chapter there is given a derivation of the leading differential equations, one or more analytical solutions, a numerical method of solution, and some results together with a dis-cussion about the applicability of the model.

2 Mathematical formulation

2.1 General Equations

In this study about mathematical models for linear simple harmonic water wave propagation problems it will be useful to have an overview of all the restric-tions that have been made in building these models. Therefore this Chapter starts with a short derivation of the basic equations of continuity and motion in Eule-rian form for an ideal (no viscosity) fluid, making use of the well-known diver-gence theorem of Gauss. The notation will be in vector form with the nabla and delta operator. '

Let V be a fixed volume in the cartesian space spanned by the three orthogonal X, y and z axes, 0 the surface area of the volume V and n the unit normal vector of 0, positive in an outward direction (see Figure 2.1). The motion of the fluid will be characterised by the scalar density p and by a velocity vector v =

(u, v, w) in which u, v and w are the components in the x, y and z direction re-spectively.

Figure 2.1 Definition Sketch

'' When £ IS a vector with the components a , a^ and a and b is a scalar, then the following vector operations will be used;

gradient of b : Vb Laplacian of b : Ab = A^b = divergence of a: V.a rotation of a : V x a = (-,8b 9b 9b, ^3x ' 9y ' 9z'' 9^b 9^b S^b 1 9x 9a, 8y^ 3a, •h-z.' 3z 3a, 3a, 3y 3z 3z 3x 3x

3y

The equation of continuity is based on the principle of conservation of mass. The change of mass M in the fixed volume V per unit time equals the change of mass in V due to the change in density p together with the outflow of mass per unit time across the surface 0.

dM _ f

dt Jj

If there is

f-[ | £ dV + f p v . n do

no production or loss of mass then

or with the divergence theorem of Gauss (/// V . a_ dV = 11 ^ • J. *iO)

V . p v) dV = 0,

I . 'If

and then, because of the arbitrariness of the volume V

^ + V . p v = 0 (2.1.1)

O t

which is the well known equation of continuity.

The equations of motion are based on the principle of Newton, laying down that the change of momentum in V per unit time equals the total force acting on V:

ƒƒ( ^ ^ ^ ^ ƒƒ P (v.il)XdO= f|[ P F d V . If K

with F representing the force vector acting on the mass in V and K the force vec-tor acting on the surface 0. When viscosity is neglected,the force acting on the surface can be represented by a pressure function p :

K = - p ri Thus 3pv p n do - \(( p F dV = 0. dV + p (v . n) v do H V "' '^•'o •'•'0 "'V

Forming the inner product with an arbitrary constant vector a^ and making use of the divergence theorem of Gauss leads to the equation

f [-""^ + V . {p (a . v) v} + V . pa - a . p F J d V = 0.

After some vector manipulation, taking into account the constantness of the vec-tor ji, the following expression can be derived.

iff a . [ - ^ + (V . p v) V + p (v . V) v + Vp - pF ] dV = 0 .

The volume V and the vector a are arbitrarily chosen so

- ^ + (V . p v) v + p (v . V) V + Vp - pF = 0.

With the aid of the continuity equation (2.1.1) and the vector rule

(v . V) V = ^ V (v . v) - V X (V X v ) ,

the equations of motion become

1= + I V (v . v) - V X (V X v) + - Vp - F = 0^. (2.1.2)

In the case of an incompressible and homogeneous fluid the density p is a con-stant, giving the equation of continuity;

V . V = 0. (2.1.3)

When the force vector F represents a conservative field,then F can be written

as

F = - Vf2,

in which ^ is the potential function of the force field. The vector equation of motion is then:

r

2•2 Energy Flux

In the following Chapters of this study a consideration about the energy flux vector will be required and therefore a general expression of the energy flux vector will now be derived.

The energy E enclosed in the volume V is given by the sum of the kinetic and the potential energy of the fluid in V:

E = [[[ [ J p (v . v) + p[j] dV.

'V

The energy flux dE/dt in this fixed volume V is then

dE

d

! = fff^^[^ p ^^-^^ ^ p " ] " '

and, in the case f^, i s c o n s t a n t in time:

Èl= [{{

A r i

d t J J ) „ 3

^ [ i p (v . v) ] dV.

Forming in Equation (2.1.4) the inner product with the velocity vector \7, and taking into account the continuity Equation (2.1.1), the following expression for the change of the kinetic energy can be derived:

^ [ J p (v . v) ] = - V . { 1 p (v . v) v} - V . Vp - pv . Vfi.

With this result and with p = constant, so V . v = 0, the energy flux becomes:

dE

jf = - III [ V . { 1 p (v . v) v} + V . pv + V . pn v j ,

and with the divergence theorem of Gauss

dE

n . n p (v . v) + p + pf21 V do.

The energy flux vector P^ per unit surface area is now defined a s :

2.3 Potential Formulation

For a potential formulation it is necessary to assume that the motion of fluid is irrotational:

V X V = 0. (2.3.1)

This assumption can be made plausible by saying that when the motion starts from a situation of rest and the fluid is ideal,there is no mechanism to introduce ro-tation. Now a velocity potential function $ can be defined from which the velo-city field can be obtained by taking the gradient:

V = V*. (2.3.2)

The field equation for the potential function $ can be obtained by the substitution of (2.3.2) into the continuity Equation (2.1.3) for an incompressible and homo-geneous fluid:

V<1> = 0 (2.3.3)

which is the well known equation of Laplace. The vector equation of motion (2.1.4) now becomes

V I?- + V {^ (V* . V$)} + V £ + Vn = 0, (2.3.4)

d L p

which can be integrated with respect to the space coordinates, yielding the equa-tion of Bernoulli:

| i + 1 ( W . V$) + £ + f2 = f(t) (2.3.5)

with an integration "constant" f(t).

The velocity field will not be affected if an arbitrary function of time is added to the potential $, so with no loss of generality the function f(t) can be chosen identically zero. In a gravity field with an acceleration g in negative z-direction,the external force potential ^ is given by

r

and then the equation of Bernoulli,

1 ^ + I (V* . v$) + £ + gz = 0 ,

(2.3.7)

frem which the pressure p can be computed.

The energy flux vector P becomes in the case of a potential flow with the aid of Equation (2.3.7)

3$

(2.3.8)

2.4 Boundary Conditions

For solving Equations (2.3.3) and (2.3.7) for the unknown three-dimensional potential function $ and pressure p,there must be given a set of conditions at the boundaries of the solution domain.

If it is assumed that the free surface elevation is given by the relation z = ri (x,y,t) (see Figure 2.2), then there must hold for a particle belonging to the free surface the following kinematic equation;

dz^ _ 9ri _^ 9n dx 9ri dy dt 3t 8x dt 9y dt

at dx 3n 3y

and with the potential function $ :

Sn _, 3$ 3ri ^ 3$ 3n 3$ 3t 3x 3x 3y 3y 3z 2 = r| (x.y.t) 0 at z = n (x,y,t) . (2.4.1)

//

^^^^^^^'^^^^^^^^^^^^.V^.r^^^^^^^^^^ .^^^^^-^^^k^v^vv^^^^^^^ z =.h(x,y)

A second condition at the unknown free surface z = n (x,y,t) is that of atmospheric pressure, which can be taken as zero with no loss of generality. So the dynamic condition at the free surface is

p = 0

or,with the aid of Bernoulli's Equation (2.3.7),

1 ^ + i (V$ . V4>) + gn = 0 at z = n (x,y,t). (2.4.2)

d t

In the case of an impermeable bottom, given by the relation z = - h (x,y), the normal velocity must be zero, leading to the equation:

3h 9h

U T T — + V ^ 7 — + w = 0 ,

dx dy

or with the potential function $:

•;^V~"*"^^":=^'^-^r-=0 a t z = - h ( x , y ) . (2.4.3) 9x 9x 9y dy 9z '-^

At other solid boundaries than the bottom also the normal velocity must be prescribed,and in the case of an unbounded solution domain in the horizontal plane some additional condition at infinity must be imposed.

Usually in diffraction problems the form of an incoming wave is prescribed and a scattered wave will be sought, which satisfies the so-called radiation condi-tion stating that the scattered wave is an outgoing and radiating wave.

2.5 Linear Equations

The non-linear boundary conditions (2.4.1) and (2.4.2), and the fact that these conditions are imposed upon the unknown free surface z = n (x,y,t),make it very difficult to find an applicable way of solving this three-simensional

problem. Therefore these boundary conditions will be linearised by means of a series development with respect to a small parameter.

There are three characteristic lengths in this problem (see Figure 2 . 1 ) :

1 : the characteristic horizontal distance of the free surface disturbance, o

D : the water depth, and

tions is

Y = H/l (wave steepness) ,

but in the case of shallow water also the parameter

a = H/D

must be small.

A derivation of the infinitesimal-wave approximation can be found in Lamb (1932), Stoker (1957) and Wehausen (1960).

The resulting equations are:

A* = 0 ( 2 . 5 . 1 ) a t z = 0 ( 2 . 5 . 2 ) 3ri_ _ 3$ 3 t 3 z 1 3 * n + - 5 - = 0 a t z = 0 ( 2 . 5 . 3 ) g 3 t and 3h ^ 3$ 3h )y 3y

5 ^ 3 l ^ ^ 3 7 f 7 ' f i = ° a t z = - h ( x , y ) . (2.5.4)

"Equations ( 2 . 5 . 2 ) and ( 2 . 5 . 3 ) can be combined by e l i m i n a t i n g the free surface

function n to the equation

| i . i ^ = 0 a t z = 0 . ( 2 . 5 . 5 )

3z g 2

2.6 Harmonic Solutions

Now that the equations have been linearised it is possible to seek a solution of Equations (2.5.1)-(2.5.4) which is simple harmonic in time;

$ (x,y,z,t) = 0 (x,y,z) cos u)t + 0„ (x,y,z) sin ujt, (2.6.1)

in which oj is the angular frequency of the wave. The simple harmonic solution (2.6.1) can also be written with the aid of a complex potential function

Re (0e ) with Re( ) indicating the real part of the ex-pression between the brackets and i = - 1 . The complex potential function 0 must satisfy Equations (2.5.1), (2.5.4) and (2.5.5) which turns into:

30

_

u^

3z B at z = 0- (2.6.3)

Once the solution 0 has been found, all the other functions of interest can be computed. The free surface function can be found with the aid of Equation (2.5.3) giving:

ri (x,y,t) = Re (— 0e )

z=0' (2.6.4)

The wave height, which is twice the amplitude of the free surface elevation, is then :

"1 -^z'O

The phase S of the free surface elevation is given by

(2.6.5)

S = - arctan (0,/0„) „,

1 2 Z=U (2.6.6)

and the linearised pressure by

p (x,y,z,t) = - pgz + up (0 sin ut - 0 cos tot) (2.6.7)

It may be useful here to summarise all restrictions of the mathematical formulation up till now. They are:

- An ideal fluid, which means no viscosity; an incompressible and homogeneous fluid; a gravity force field;

an irrotational motion, which leads to a potential formulation; - infinitesimally small amplitude waves; and

r

2.7 Fundamental Solution

Assuming the bottom is horizontal, the bottom condition is

£*i = 0 at z = -h. (2.7.1)

dz

A solution of the Laplace Equation (2.5.1) with the boundary conditions

(2.6.3) and (2.7.1) in the z-direction can now be found by the method of separation

of variables. By writing

0 (x,y,z) = Z(z) <}. (x,y), (2.7.2)

and substitution into the Laplace equation, it can be derived that

1 i ! i = - 1 (i!i + l ! i

(^^^ + ^ ^ ^ ) , (2.7.3)

dz^ '*' 8x2 gy2

giving the two equations

} ^ = b 2 (2.7.4)

dz

and

* 8x2 gy2

where b is a complex constant (eigenvalue).

The solution of Equation (2.7.4), satisfying the free surface condition and the

condition at the bottom, is given by:

Z (z;b) = C(b) cosh {b (z+h)}, (2.7.6)

provided that the eigenvalue b is a root of the dispersion relation

u^ = gb tanh (bh) (2.7.7)

which can be derived by substitution of (2.7.6) into the boundary condition (2.6.3).

Equation (2.7.7) has the real roots b = ± k and an infinite number of pure

imaginary roots b = ± ik.; j = 1,2,3, (see Figure 2.3).

2TT 3TT _ ^ k h - t a n h (.kh) u'h _L 9 kh 2 3

— * kh

Figure 2.3 Roots of the Equation ti)^ = gb tanh (bh)

The solution can now be written as a superposition of all eigensolutions, giving (x,y,z) = Z Z (z) (t).(x,y) j=0 J J (2.7.8) with Z.(z) = C. cosh {T.(z+h)} (2.7.9) T = k ; T . = i k . , j = l ,2,3, o o J J -^

and the two-dimensional potential functions (J).(x,y) which are solutions of the equation

— 1 + J + T.^ <t>. = 0. J J

(2.7.10)

The eigenfunctions Z.(z) form a set of orthogonal functions in the interval - h ^ z-^ 0, which can be proved with the aid of the field Equation (2.7.4) and the boundary conditions at the bottom and free surface. From Equation (2.7.4) it can be shown that

(T.^

-

T

^) r z.

J m J J

-h Z dz

r°

d^z. d^z

(Z i - Z. 5i) dz,

(

"> dz^

J

dz^ -h

r

(T.^ - T 2) J m r (- dZ. Z. Z dz = Z ^-J- - Z. J J m [_ m dz j dZ 1 Tl dz z=0 z=-h

and with the boundary conditions

dZ. dZ 2

^-J- = 0 at z = -h and - ^ = — Z. at z = O dz dz g j

for all j = 0,1,2,.... there can be found O

(T.^ - T ^) [ Z. Z dz = 0, J m J J m

-h

giving the orthogonality property 0

I Z, Z dz = 0 for ra ?* i.

J J m (2.7.11)

The coefficients C. can be chosen arbitrarily because the function *. in the

J J , solution (2.7.8) can take up also constant factors. These constants C. are taken in such a way that the eigenfunctions Z.(z) now form an orthonormal set of functions in the interval - h ^ z ^ 0.

0

Imposing Z. dz = 1

gives the result

C.^

h (T.-" - V) + V

for j = 0,1,2, (2.7.12)

in which v = to /g.

2.8 Solution with Singularity Distributions

A solution of the three-dimensional problem of wave propagation in a certain domain with obstacles and a variable water depth can often be found with the aid of a singular source solution G(P;M), which possesses the following properties:

- As a function of the coordinates of the point P the source potential

G(P;M) is a solution of the Laplace equation except at the source point M;

the source potential satisfies the free surface condition and the flat

. . 30

bottom condition 77— = 0 at z = -h;

dz

- G(P,M) satisfies the radiation condition at infinity, which is given by

lim

JT

(1^ - ik G) = 0 (2.8.1)

r-j-to

with r the horizontal distance between the points P and M;

- G(P,M) is defined and is regular analytic in the whole domain except at

the source point M; and

the singularity in G is of the form 1/R with R the distance between the

two points P and M, so the function G - 1/R is regular analytic in the

whole domain.

Looking at the fundamental solution given in section 2.7, an expression can

be found for the source potential function G(P;M), which is also symmetrical with

respect to the coordinates of P and M:

G(P;M) =

T.

C.2 cosh {T.(z+h)} cosh {T.(C+h)}

4^ H

^(T.r) (2.8.2)

j=0 J J J 2i o J

with (x,y,z) = the coordinates of the point P

(?.r|jC) = the coordinates of the source M

r =

s/

(x-C)2 + (y-n)"'

H (Tr)= a solution of the two-dimensional Equation (2.7.10), satisfying

the radiation condition (2.8.1) and possessing the apropriate

sin-gularity for r = 0 (Hankel function of the first kind and zero'th

order),

By separating the first term in the series development of (2.8.2), the source

potential G can be written in a somewhat different form, taking into account the

expressions x. = i k. ; i = 1,2,3,.... and T = k :

J J o o

2(k^-v2)

G(P;M) = cosh {k (z+h)} cosh {k (C+h)} ^r^ H ' (k r) +

o o 2i o o

, 2(k.^+v2)

- - E

i

cos {k.(z+h)} cos {k.(C+h)} K (k.r) (2.8.3)

^ j=l h(k.^fV^)-v J J o J

with K (k.r) = Modified Bessel function of the first kind and zero'th order

There exists also an integral expression for the function G(P;M), showing the singularity like 1/R at the source point M:

2TTG(P;M) = - - - — + ^ R' cosh (kh) {k tanh (kh) -v} o

2ir(k 2-v^)

i cosh {k (z+h)} cosh {k (C+h)} J (k r ) ( 2 . 8 . 4 )

1 /, 2 2\ O O O O h(k ^-v )+V with R = y (x-C)' + (y-n)' + (z-C)'' ^ / ( x - C ) ' + (y-n)' + (z+<;+2h)2

R'

J ( ) = Bessel function of the first kind and zero'th order o

P.V. = The principal value of the integral, because of the singularity of the integrand for k = k .

The first term in the right-hand side of Equation (2.8.4) is a general solution of the Laplace equation for an unbounded domain. The second term is needed to satisfy the condition at the flat bottom and the other terms are required for fulfilling the free surface condition at z = 0.

More about singular solutions and their representations by series developments and integral expressions can be found in articles by John (1950), Thorne (1953) and Wehausen (1960).

The diffraction problem in three dimensions, which is to find the scattering of an incoming wave due to the presence of obstacles and a variable water depth, can be solved theoretically by distributing wave sources over the surfaces of the ob-stacles and variable bottom which disturb the incoming wave (see Figure 2.4):

0(P) = 0(?) + j a(M) G(P;M) dOj^ (2.8.5)

0

in which 0 = the known wave potential of the incoming wave, a( ) = the strength of the source distribution.

z = ^ ( x . y . t )

, P ( x . y , z )

V

H\Wy.^-JM^V/f^

W/'AVAWAWAV/AV-'AWAV)

Figure 2.4 Three-dimensional Source Distribution

The strength a of the source distribution can be found by the condition of re-flection at the disturbing surfaces, giving the integral equation:

a(P) +

a(M)

3G (P;M)

3n •^"M = - ^3^>P> (2.8.6)

with both points P and M situated on the surface.

When a solution of this integral equation of the Fredholm type has been obtain-ed the desirobtain-ed solution 0(P) at any point P can be computobtain-ed with expression

(2.8.5).

This method of solution for wave propagation problems has been used by Garrison (1972) and v. Oortmersen (1973) to compute the wave forces upon submerged obsta-cles due to an incoming wave. They solved the integral Equation (2.8.6) numeri-cally using a technique that is very similar to the panel method in aerodynamics; Hess (1971). However, to use this method for computing the wave propagation over a large area, the computing time and computer storage will be excessive and there-fore there is need for some simpler mathematical models restricted to the horizon-tal plane.

3 Diffraction

3.1 The Diffraction Equation

A reduction to a two-dimensional diffraction model can be obtained when the variation of the bottom profile can be neglected (the mean slope of the bottom over a certain characteristic distance, for instance»the water depth, is zero) and when the obstacles are cylindrical extending from the bottom to the free surface.

In that case the solution can be written as a superposition of a known incident wave potential 0 and a scattered wave potential 0 which has the form (2.7.8) as has been derived in Section 2.7. The two-dimensional Equation (2.7.10), the

re-30^ 30

flection condition -;r— = - -jr— at the contours of the obstacles, and the radiation 3n dn

condition define the functions (J).(x,y) ; j = 0,1,2

The incident wave is given as a plane wave and has the form: -? cosh {k (z+h)} ^

0(.,y,.) = - I - ^ i ^^^^l^ ^(x.y), (3.1.1)

o ^ i k (x cosg + y sing)

in which ((i(x,y)=e (3.1.2)

corresponding to a free surface elevation

n(x,y,t) = I H cos {k (x cosg + y sing) - a)t} , (3.1.3)

With H = wave height of the incident wave.

3 = the angle between the direction of propagation and the positive x-axis. k = the wave number belonging to the water depth h and the angular

fre-quency O), accordingly to the dispersion relation (2.7.7).

The reflection condition — — = - -?T— at the vertical walls can only be fulfilled 9n dn

for all values of z, when

" 5 ^

^ - 0 for j = 1 ,2, (3.1 .4) and

at the contours of the obstacles.

It can be shown now that the functions (j>.(x,y) for j = 1,2,.... must be identi-cally zero. From Equation (2.7.10) it can be derived that for the area A (see Figure 3.1) the following equality must hold (the nabla operator is now two-dimensional) :

r ( V . ((). V(t>.) dA - r r [(?<(). .^<i,.) + k . ^ ip^idA = o .

or, with the divergence theorem of Gauss,

•'j 311

r+c

^ d s - I [(V*. .V*.) . k . ^ 0.^]

dA = 0.

Figure 3.1 Diffraction in the Horizontal Plane.

Taking now the limit r-'•'^ and taking into account that i^ . is an outdying radia-ting wave, which also satisfies condition (3.1.4) at F, there remains;

ff [m..V(i).) + k.^ (ti.n dA = 0.

Because of the fact that the integrant is positively definite, there must hold

.(x,y) = 0. (3.1.6)

poten-tial c|) needs to be computed. The differenpoten-tial equation for this function is thus:

3!i + 3!i + k 2 <j, = 0, (3.1.7)

which is the so called Helmholtz equation.

3.2 Boundary Conditions

It has been seen in the foregoing Section that the boundary conditions for the two-dimensional scattered wave potential are:

~^= -|i (3.2.1)

at full reflecting contours F and a radiation condition

r- ^*d

lim V^r

i-j-

- i ^o 'J'd^ " ° (3.2.2)

r->- C O

at infinity. With these conditions the scattered wave potential will be uniquely defined.

In practical problems the walls of the obstacles are not always vertical and often also not fully reflecting. To treat these boundaries in this diffraction model the walls must be schematized as vertical, but it is possible to introduce partial reflection into the model by using a mixed boundary condition for the to-tal wave potential c|) = {[) + Ó, instead of the full reflection condition TT— = 0.

d dn Assume the condition

^ + k a (f) = 0 (3.2.3) on o

at the partial reflecting boundaries, in which a (= a +ia„) is a theoretical non-dimensional complex reflection coefficient.

With the aid of the expression (2.3.8) for the energy flux vector P it can be shown that with this mixed boundary condition a mean normal energy flux through the boundary will be imposed. The mean energy flux P, defined by

in which T is the period of the wave, is expressed in the potential function 0:

P = ^ p 10 i (0 V0 - 0 V0) (3.2.4)

with an overbar indicating the conjugate complex value. Use of the expression , ^ cosh {k (z+h)}

0(x,y,z) = - 2 i f ^ cosh ( k h ) - '^('^'y)

O

for the two-dimensional diffraction model and computing the integrated mean ener-gy flux P , defined by

0

ƒ

-h

gives the expression

P*=Yg- P 1 ^ Ujj i (()) W - * V * ) , (3.2.5)

o

with

" E = 2 " ^' * sinh (2°k^h)^ ' S " " P velocity

c = (jj/k ; phase velocity

and the nabla operator V now two-dimensional. The mean integrated energy flux normal to the boundary becomes, with the aid of the mixed boundary condition

(3.2.3),

- * 1 ^ 2

<^B • Z ) = - Yg- P Hg U ^ i (a - a ) (J)(i) ,

or with i (a - a) = 2a and ^^ = H^/H , accordingly to definition (2.6.5) of the wave height H,

(n . P *) = - -^ p g U^ a^ H 2 . (3.2.6)

When the imaginary part of the theoretical reflection coefficient a is positive, there is an energy flux at the boundary into the negative normal direction. This absorbed energy flux is relative to the square of the wave height, assuming the reflection coefficient a is independent of the potential function (|>.

The mixed boundary condition (3.2.3) can also be translated into two conditions for the wave height H and the phase function S, which are defined by the expression

<f) = ^ e^ ^ (3.2.7)

H

These conditions are:

- 4 1 ^ = ^ a, and - 1 ^ = k a„. (3.2.8)

H dn o 1 dn o 2

It can be seen that when a- ?^ 0 the lines of equal wave height make a non-zero

angle with the normal at the boundary and that when a„ ^ 0 the lines of equal

phase make a non-zero angle with this normal.

In practical hydraulic work the reflection coefficient is introduced in a more

physical way and is defined as the ratio between a reflected wave and an incident

wave. This definition is only useful in the case of reflection of long-crested

plane wave against an infinitely long breakwater, assuming partial reflection has

no influence upon the plane wave form of the reflected wave. In that case it is

possible to obtain a relation between the theoretical reflection coefficient a

1 R

(= a +ia ) and the practical coefficient p (= Re ) in which

R = amplitude reduction factor, and

B = phase shift.

Assume a one-dimensional wave is partially reflected against a vertical wall

at X = 0. The total potential of the wave field is written as the superposition

of the incoming and a reflected wave;

(i> = j n e

° + ^ H R e ° (3.2.9)

At the reflecting wall the theoretical condition is:

TT- + k ai> = 0 a t x = 0 ,

dx o

and substitution of the expression (3.2.9) into this condition gives the two

equations:

a (1 + R cos g) - a R sin g = R sin B,

and

a R sin B + a (1 + R cos B) = 1 - R cos B,

In Figure 3.2 the values of a. and a^ are given as functions of the amplitude reduction factor R for different values of the phase shift

3-Figure 3.2 Relation between Theoretical and Practical Reflection Coefficient

The mixed boundary condition is only applicable in the case that the theoreti-cal reflection coefficient a is not dependent on the solution (J) (linear condition). However, in general partial reflection is a non-linear phenomenon and therefore the linear formulation can be wrong in the neighbourhood of partially reflecting walls. But it might be possible that with this formulation of energy absorption reasonable overall results can be obtained when the interest is only for an aver-aged wave height over some area.

3.3 Some Analytical Solutions

To get some idea of the diffraction phenomenon two analytical solutions will now be given, one illustrating the scattering of an incident plane wave field by a circular pile and the other showing the diffraction by a semi-infinitely long fully reflecting vertical breakwater.

been given by McCamy and Fuchs (1954) with the aid of a Bessel series development.

t ;.PCre.z)

Figure 3.3 Definition Sketch of Circular Pile

If (r,0,z) are the cylindrical coordinates of the point P, then the solution in this point is given by:

c o s h {k ( z + h ) }

5 ( r , 9 , z ) = - i 2 g

2 0) c o s h (k h) (t>(r,9) w i t h ( 3 . 3 . 1 ) -K r , e ) = S e i A (k r ) c o s (mO) „ m m o m=0 ( 3 . 3 . 2 ) i n w h i c h e = 1 ; e = 2 f o r m = 1 , 2 , 3 , A (k r ) = m o J ( k r ) Y ' (k a ) + J ' (k a ) Y ( k r ) m o m o m o m o J ' (k a ) + i Y ' (k a ) m o m o

J ( ) = Bessel function of the first kind and m-th order,

m ' Y ( ) = Bessel function of the second kind and m-th order,

m ' a = radius of the pile.

and a prime indicating a differentiation with respect to the argument. The wave height pattern and the phase lines are shown in Figure 3.4a and Figure 3,4b for the case k a = 1,

Figure 3.4a Diffraction around Cir-cular Pile; Lines of Equal Wave Height every 0.1 unit.

Figure 3.4b Diffraction around Cir-cular Pile; Lines of Equal Phase every

7T/4 rad.

In front of the cylinder a reflection pattern can clearly be seen and behind the pile the influence of the diffraction effect.

The problem of diffraction around a semi-infinite breakwater has been solved by Sommerfeld and a numerical evaluation of the solution can be found in Shou-Shan Fan (1967) and Born (1965). The incident wave is making an angle of B degrees with the positive x-axis (see Figure 3.5).

y y y / • y y ^ incident w

/

/

e /

y\ \ breakwater » X

Figure 3.5 Definition Sketch of Semi-infinite Breakwater

!Z;(r,e,z) = •S cosh {k (z+h)}

1

Hg 2 O) cosh (k h) <l>(r,e). (3.3.3) with (i>(r,e) = f(aj)e i k r cos(9-6) i k r cos(e+6) in which 0 = 2 Hö^)i (3.3.4) k r TT sin j (9-3) k r a, = - 2 ,/ - ^ sin 1 (8+6)

and f (a) = i ^ f ^-^^'i'l^ dq.

The wave height pattern in the case of a normally (g = 270°) incident wave is sho\m in Figure 3.6 (Shou-Shan Fan) and more in detail in Figures 3.7 and 3.8 (Born) There can be distinguished the region of reflection in front of the breakwater, the region of shadow behind the breakwater, and the region of transmission, where wave heights greater than one can occur due to diffraction effects. The phase function shows undefined points in the reflection area corresponding to points with a minimum in the wave height. An explanation of the existence of these points can be found in a report of Radder (1974).

Figure 3.6 Diffraction by a Semi-infinite Breakwater; Lines of Equal Amplitude (after Shou-Shan Fan)

Figure 3.7 Diffraction by a Semi-infinite Breakwater; Lines of Equal Amplitude (after Born)

Figure 3.8 Diffraction by a Semi-infinite Breakwater; Lines of Equal Phase (after Born)

3.4 Wave Penetration into Harbours

The mathematical model for computing the wave penetration into a harbour of arbi-trary shape is based on finding the solution of the diffraction equation (3.1.7) in an area which is bounded by vertical barriers and an aperture to the open sea (see Figure 3.9).

Figure 3.9 Definition Sketch of Harbour Diffraction Problem

A reflection condition can be imposed at the barriers, while at the entrance of the harbour the incident wave is given and the radiation condition also taken into account.

For the method of solution it is convenient to distinguish two areas: Area I inside the harbour and Area II outside the harbour. Both areas are assumed to have the same constant water depth and therefore also the same wave number k . o The boundary between the two areas is indicated by C. In Area II the

solution will be written as a superposition of the known incident wave (fi and a diffracted wave (j), caused by the presence of the harbour and coming only from the

aperture C. In this outside area reflections from the breakwaters and the coast are ignored, so the total potential function in this area is:

* = ^ + 4)^ . (3.4.1)

Inside the harbour the solution will be indicated by the potential i);. Both (J) andii will be sought with the aid of singularity distributions. The source solution of the Helmholtz equation satisfying the radiation condition is given by

G(P;M) = ^ H '(k r) (3.4.2) zi o o

The diffracted wave (fi in Area II is written as a superposition of source potentials coming from points situated only at the aperture C.

<1>,(P') = r y^/M') -!v H ^ k r') ds (3.4.3)

d _/ II 2i o o f^i ^,

C

The inside potential ijj in the point P can be written as:

i^m =

ƒ Pj(M) 2jH^'(k^r) dsjj. (3.4.4)

r+c

In these expressions y (M') is the intensity of the source M holding for Area I and y (M) is the intensity of the source M holding for Area II (see Figure 3.9) Imposing the reflection conditions at the boundary F and requiring continuity conditions for the normal velocity and wave height at the aperture C, the unknown intensity functions y and y must satisfy the following integral equations:

(|^)_ = ^ i ^ ^ ) ^ { / i 3 ! ^ [ ï r V ( V ' ] < ^ ^

p c+r

and

c+r c

for P situated on the aperture C and

I J I dn L zi o o J J I 2i o o ei A T'i

C+r r ( . j - ^ - / ^

for P situated on the reflecting contour F with the theoretical reflection coef-ficient a (see section 3.2). With these integral equations the source intensity functions \i and \i are uniquely defined and can be solved numerically.

Once the intensity functions have been found the potential at each point in Area I or Area II can be computed according to the expressions (3.4,4) and

(3.4.3).

For complicated harbours with many basins (see Figure 3.10) it is possible to split up the harbour into more areas corresponding to the basins and to express the solution in each area as a source integral over the boundaries of that area.

All—'

Figure 3.10 Harbour with Multiple Connected Basins

Requirements of continuity for the normal velocity and wave height at the bound-ary between two areas create a set of integral equations for the unknown inten-sity functions y of all basins.

3.5 Harbour Oscillation

The mathematical model for wave penetration into harbours of arbitrary shape can also be used for the computation of harbour oscillation. In that case the

interest is in the wave heights at several points of the harbour as a function of the period of a disturbance at the entrance. For some periods the wave height shows a maximum and the corresponding frequencies are called the resonance fre-quencies of the harbour. These resonance frefre-quencies are not much effected by the form of the disturbance at the entrance,but the value of the wave height depends on the incoming normal energy flux and therefore on the width of the entrance. As an illustration of the harbour resonance phenomenon,a simple example of the oscil-lation modes of a long rectangular harbour with a full open entrance is given in Figure 3.11.

r

It is assumed that the wave form in the harbour is nearly one-dimensional, so in

Area I the wave potential is given by:

i k X -i k X

(t)^(x) = Ae ° + Be ° , (3.5.1)

with A and B integration constants.

In the area outside the harbour the wave potential (|) is a superposition of a

CO

known disturbance (() and a diffracted wave (j) , which can be written as a source

integral over the entrance C:

+ Jb

*jj = * + ƒ y(y) 21 H Q ^ V ^ ''y ' (3.5.2)

-ib

with b the width of the harbour entrance. The reflection condition - i — = 0 at the

an

end X = L of the harbour, together with the two continuity conditions:

a * 3<t>

<i>^ = ^^^ and ~ = ^ ^

(3.5.3)

^I ^11 an dn

at the entrance x = 0, give three equations for the three unknowns A, B and y.

From the reflection condition there can be derived:

(() = C cos {k (x-L)}, (3.5.4)

with C a new integration constant. The two continuity conditions (3.5.3) give the

equations

+ lb

?| + ƒ U(y) 21 "o'^'^ol^l^ "^^ " ^ ''"^ ^""o^^ (3.5.5)

-jb

and

C k sin (k L) = - -^1

+ \i .

(3.5.6)

o o dx

'x=0

Assuming the disturbance c|) is an incident wave travelling in positive

x-direc-tion and that there is no reflecx-direc-tion outside the harbour, then

and

I't 1 - i k X

| i = - l H i k e °

d x 2 o

The e q u a t i o n s f o r t h e unknown C and y become

Dp - C c o s (k L) = - TT H o I and 1 <^ u + k C s i n (k L) = - 77 i k H , o o 2 o w h i c h h a v e t h e s o l u t i o n s 1 - i k D ,•<> c o s (k L) + k D s i n (k L) 5H 0 0 o ( 3 . 8 . 7 ) s i n (k L) + i c o s (k 1)

y ^ _ o ° n 8 o-,

, , -ï? c o s (k L) + k D s i n (k L) ^ • • ) 5 k H 0 0 o o

with D = ƒ Yi H o ' ( k j y l ) dy.

The c o e f f i c i e n t D c a n b e a p p r o x i m a t e d f o r s m a l l v a l u e s of k by b r 2

D

.[1 (Y_, + In (Ik^ b)} - i] , (3.8.9)

with Y = 0.5772157 (Eulers constant).

Figure 3.11 shows the value of C/^H as a function of the wave length

for different values of the width b. In the limit b -> 0 the resonance wave

lengths X, are given by

Both the three-dimensional and the two-dimensional diffraction problem lead to an integral equation for the source intensity function over the surface or contour of the obstacle. As these integral equations cannot be solved analytically for surfaces or contours of arbitrary shape, a numerical method must be applied. For three-di-mensional problems reference can be made to the work of Garrison (1972) and for

two-dimensional problems to Daubert and Lebreton (1965).

The method will be illustrated in this section for the two-dimensional problem of diffraction around a single cylindrical full-reflecting obstacle. The method is also applicable with slight modifications to the problem of wave penetration into harbours with partial reflecting walls. The numerical method of solution can be classified as a finite element method because of the following three steps: - The contour of the obstacle will be split up into a number of finite segments;

in each segment the unknown intensity function will be approximated by a func-tion with some parameters; and

- continuity conditions between the approximation functions in the connection points together with the integral equation define uniquely the unknown para-meters of the whole problem.

\

r

^

y

t

• * — S r/

/ J

-^ - P ( a )

Figure 3.12 Definition Sketch Method of Solution of Two-dimensional Diffraction Problem

The integral equation for this problem is given by:

,(P) .ƒ,(„) _|_ [^H^(V)]^3„=-(g)

(3.6.1)

and after splitting the contour F into a finite number of segments P. ; j = 1,2, , N; by:

N

y(p) + r

ƒ "

r.

j

(M) ^ r ^ H i(kr)]ds^ = -(!*)

diL, L 2i o o J M ^on'o (3.6.2)

This integral relation will now be imposed in a finite number of discrete points P ; m = ! , 2 , ,N (collocation points), which are the centre points of the

seg-m

ments (see Figure 3.12):

V(^J

N

V(M) 3 r ' H >(k r)] ds„ = - ( | i )

9n^ L 2i o o J M ^'^"'p

(3.6.3)

Not only must the intensity function p be approximated piecewise but also the form of each contour segment to carry out the integration in expression (3.6.3). At this stage both these approximations must be attuned to each other, because it is useless to obtain an accurate solution belonging to a bad contour approximation. In general, both approximations must be of the same order. Let the form of the contour-segment r. be approximated by a polynomial expression (spline) of some order. De-pending on the order of this spline, continuity conditions for the derivatives be-tween the splines can be imposed to obtain a smooth approximation of the curve T. In the same way the unknown intensity function in segment T. will be approximated by the spline Li,(s:a. , a.,, ,a. ) of order p, in which a. , a.,, etc. are

^ J Jo' jl' ' jp _ ^' JO* jl'

the p+1 parameters of the spline. With a cubic spline (p=3) continuity conditions for the intensity function itself as well as the first and second derivates can be imposed in the connection points, by which only one unknown parameter per segment remains. These parameters can now be computed with the aid of the discrete inte-gral relation (3.6.3), which gives a set of linear equations for the unknown para-meters.

It is obvious that by taking higher order splines for the approximation of both the contour ,3nd the intensity function along the contour a lot of work must be done to compute the coefficients of the ultimate set of linear equations. In gen-eral, the integration over each segment must be done numerically because of the complicated source function. On the other hand, the gain of all this work can be a saving in the number of segments and therefore a saving in computing time to get a solution of the set of linear equations.

The simplest way to do the numerical computation is to split the contour in-to a number of straight~line segments and taking both the intensity function and the source function in expression (3.6.3) as a constant. The work of computing the coefficients of the system of linear equations is now minimal and also the approximation of the contour is very simple. Of course, sufficient segments must be

taken to approximate the curved contour m an appropriate way and to adequately represent the exact intensity function in discrete points.

In practical computations this most simple numerical method of solution will be used because for practical problems such as wave penetration into harbours the requirements for the accuracy of the results are not severe, taking into account the uncertainty of the input data of the model.

An indication of the influence of the number of computing points upon the re-sults can be obtained from some computations in the case of diffraction around a circular pile (see Figures 3.13a and 3.13b for a three-dimensional computation and Figure 3.14 for a two-dimensional computation).

Figure 3.13a Amplitude of the Wave Po- Fig. 3.13b Phase of the Wave Poten-tential 0(x,y,z) along the Surface of a tial (|)(x,y,z) along the Surface of a Pile for Different Schematizations. Pile for Different Schematizations. n = number of points along the

circu-lar contour

n = number of points along the verti-cal

It seems that a minimal number of computing points over a distance of one wave length has to be maintained to ensure reasonable results. For practical problems this number must be five or more.

About the influence of higher order spline approximations the reader is re-ferred to the work of Botta (1975) and Labruyere (1972), both working with similar

problems in the field of aerodynamics. H / H - - e x a c t • simple method + constant intensity polygon appr

O

« quadratic

intensity-A constant i n t e n s i t y , circular appr f j

:-4

> numerical 6 segments

Ê)

y

< ' " - . - 1 , 4 '

' \

Figure 3.14 Influence of the Number of Computing Points upon the Wave Height at r = 30

3.7 Some Results

It is not the purpose of this study to give very accurate and detailed solutions of some diffraction problems, but more to illustrate the possibilities of the diffraction models for practical problems.

First the results of a harbour diffraction computation are compared with some measurements in a hydraulic model, in which the walls were vertical and fully reflecting, the bottom horizontal and the incoming wave sinusoidal, to ensure the best possible simularity with the mathematical model. The lay-out of the

harbour can be seen in Figure 3.15 and also the measuring-points in the harbour. In Figure 3.16 the wave heights, relative to the incoming wave height, at the measuring-points are given for both the computation and the measurement. The agreement between measurement and computation is very good, especially the tran-sition to higher wave heights. The computed wave heights, however, can deviate

The last result belongs to a harbour resonance problem (seiche) for a harbour of arbitrary shape. In Figure 3.23 the relative wave heights at points A and B of the harbour are given as a function of the period of the incoming disturb-ance at the entrdisturb-ance. The walls of the harbour are fulfy ref1ecting and the dis-turbance is given as a plane wave normal to the entrance. In this Figure the wave heights show some peaks, and the corresponding periods are called the periods of resonance. When these periods are out of the range of periods for all possible disturbances at the entrance, then the harbour will have no resonance problem. Be-cause of the longer periods, which will be considered, the number of computing points along the contour is not very high, thus allowing a repetition of the compu-tation for more periods within a reasonable computing time.

Figure 3.15 Lay-out of the Harbour

ncident wave

f-"'"-^ /

\

Figure 3.17 Wave Height Contours

incident wave

Figure 3.19 aj = O, 82 = 1/3 Figure 3.20 aj = 0.4, a2 = 0 Figure 3.21 a = 0.4, a^ = 1/3

Diffraction Around a Partial Reflecting Circular Island Lines of Equal Wave Height every 0.1 unit

^ " ' \

Figure 3.22a Diffraction Around a Partial Reflecting Circular Island (a. = 0, a = 1 ) . Lines of Equal Wave Height every O.I unit

Figure 3.22b D i f f r a c t i o n Around a P a r t i a l Reflecting Circular Island (a. = 0,

a = 1 ) . Lines of Equal Phase every TT/4 r a d .

50

100

150

200 250 300

> period (sec)

3 5 0

4 0 0

4 5 0

500

3. 8 Discussion

In spite of the restrictions for the diffraction models and mainly the require-ment of linearity, it appears that the results give qualitatively and in some cases also quantitatively reasonable answers to the problems. Of course, the phenomenon of separation of flow and the creating of vortices are not included in this linear mathematical model»and therefore the results for structures or harbours with

sharp corners will be less accurate, especially in the neighbourhood of these points.

In the case of wave force computations Garrison (1972) made some comparisons with measurements on a smooth obstacle and found reasonable results.

In harbour diffraction problems the lines of equal wave height correspond very well to the measurements, and for practical problems of harbour design this will be sufficiently accurate. Requiring more accurate results will not be meaningfull because of the uncertainties in the boundary conditions in nature, such as period and direction of the incident wave. Also the phenomenon of partial reflection of the breakwaters gives rise to the problem of choosing the appropriate reflection coefficients. Partial reflection has been observed in the past mainly for perpen-dicular incident plane waves. Measurements show that partial reflection, in which the breaking of the waves is included, is a non-linear phenomenon, depending on the slope and roughness of the breakwater and on the wave length and wave height near the breakwater. This means that the theoretical reflection coefficient in-troduced in Section 3.2 will depend in general on the wave potential. It might be possible, however, that the linear formulation as a mechanism of energy dissipa-tion along the boundaries will give reasonable overall results,although this must be shown by extensive measurements in a hydraulic model. Also there is only a little knowledge of reflection coefficients in the case of obliquely incident waves, and the aspect of the phase shift due to reflection has had little attention by workers in this field.

Results of harbour oscillations computations have been compared by Lee and Raichlen (1971) with measurements in a hydraulic model and the agreement in reso-nance periods is good. Of course,the values of the wave height belonging to these resonance frequencies are not realistic, because loss of energy in the so-lution domain is not possible in this mathematical model.

4 Refraction

4.1 The Refraction Equations

The phenomenon which is called as refraction of waves can be described as the change in the direction of wave propagation, in the wave length and in the wave height due only to a slow variation of the bottom. Reflections of the waves against boundaries in the horizontal plane can be ignored.

The differential equations which describe the refraction phenomenon can be derived from the basic linear potential equation (2.5.1) and the conditions at the free surface (2.5.5) and at the bottom (2.5.4) by means of an asymptotic ap-proximation. In this theory it is assumed that, corresponding to a slow variation of the water depth,the variation of the potential in the horizontal plane is much less than that in the vertical direction. The refraction equations de-rived with this theory are essentially the same as the classical equations, which can be derived in a geometrical way using the laws of Snellius. The asymptotic approach used in this Section differs slightly from that of Keller (1958), who assumes in advance that the phase function is only a function of the horizontal coordinates.

Assuming a slow variation of the bottom, the gradient of the water depth in the bottom condition (2.5.4) can be made of the order one by taking new horizontal co-ordinates (x,y) in the water depth function h(x,y):

(x,y) = ^ (x,y)

in which (see Figure 1.1) D = the mean water depth, and

a = the mean slope over a distance D,

giving:

Vh = ^ Vh (4.1.1)

with V and V the two-dimensional nabla operators (7—- , -;--) and (•;?= , 77=) respect-9x 3y respect-9x 3y ively.

The other coordinates are made dimensionless with the aid of the free surface characteristic length 1 = g/u^

(x', y', z'.h') = y- (x,y,z,h). o

The equations written in these quantities, omitting the primes for simplicity in notation, are now:

i!l.Ü^.Ü^=0, (4.1.2)

8x2 3y2 32^

1^ - 0 = 0 at z = 0, (4.1.3)

dz

8 É . , ( | « | h ^ 3 0 3h ^ ^ _

3z dx 3x 3y 3y

with e = a/u and y = D/1 .

o

Essential for the refraction approximation is the assumption of slow variation

of the potential in the horizontal plane. Mathematically the horizontal

coordi-nates X and y in the potential function are scaled with the parameter e, giving

the equations:

,2 ( Ü 0 . Ü 1 ) , 3 ! 0 . O , (4.1.5)

9x2 3y2 3z2

1^ - 0 = 0 at z = 0, (4.1.6)

dz 30 2 ,30 9h ^ 30 3h. _ - ,^ , - V, ri , 7-, ^ ^ ^ ^3Ï 35 ^ ^ 3 f ^ - ° at z - -h, (4.1.7)

in which it is assumed that all quantities, except the parameter G, are of the order one.

There is now introduced an amplitude function A(x,y,z) and a phase function S(x,y,z) by writing:

0 = A e"- ^.

-2 ri

z^

{ | V^A - (VS .VS)} ,

e^ V.A^VS + ^ (A^ 1^) = O,

dz dz 1 8^A ,8S

(ff)^ = °'

(4.1.8) (4.1.9)

dz

at z = O, (4.1.10)

1^ + E^ (VA . Vh) = 0

dz (4.1.11)

3z

at z = O, (4.1.12)

3z

e^ (VS . Vh) = 0 (4.1.13)

Equation (4.1.9) is the equation for the energy flux and can be integrated over the water depth,thus:

V . A^VS dz + A^

3z

z-0

z=-h

or, using the boundary conditions (4.1.12) and (4.1.13),

V . A^VS dz + A^ (VS . Vh)| = 0 -h 'z=-h

ƒ

V . A'^VS dz = 0. -h (4.1 .14)

Now solutions of the functions A and S will be sought in the form of series de-velopments with respect to powers of the parameter e. An approach with the aid of order functions, as is treated by v. Dyke (1964), is more general, but this method will complicate the derivation and will give no other results. It appears that the series development is in even powers of G, so:

A = A + e^A, + e'A^

and

S = E 1 2

fs + c^S, + e''s„

L o 1 2

(4.1.15)

with m a power, which must be defined later on. Substitution of these develop-ments into Equation (4.1.9) and condition (4.1.12) gives, in first approximation

for small values of e, ~ 3S dz O dz

or,after integration and using (4.1.9),

^ = 0 , (4.1.16)

which means S = S (x,y). Using these results and substituting the series devel-opments in equation (4.1.8) gives:

e^V^A - e^^*^ÏA (VS . VS ) + e^ {A, (VS . VS ) + 2A (VS . V S , ) } 1 +

_o

_L

_o

_o

_o

_l

_o

_o

_o

_o

_l

_j

3^A 3^A, ^ 3S,

+ e^ e A (~—)^ + 0 ( e ) + 0(e ) = 0.

3z2 9z2

2 O 9z

Depending on the value of the power m,three possiblities can be distinguished. First,m > O gives in first approximation the equation

3^A

—

- =

0,

and with the aid of condition (4.1.11)

9A

or A is only a function of the horizontal coordinates x and y. However, condi-o

tion (4.1.10) requires then the only possible trivial solution A = 0. In the second case, m < - 1 the first approximation is:

A (VS . VS ) = 0, 3S c

"3x

function as the only possible trivial solution.

The remaining possibility is -1 < m < 0, giving the equation: 3^A - A (VS . VS ) = 0 ^ 2 O O O dz , 3^A — = (VS . VS ) . o 9z

As the right-hand side of this equation is only a function of x and y, so also must the left-hand side be a function of only x and y. Let

(VS . VS ) = K^(x,y), (A.) .17)

o o

then the amplitude function A must be a solution of the equation

a^A

- K^A = 0 (4.1.18) 3z^

with the boundary conditions

3A ^ = 0 at z = -h (4.1.19) dz and 3A - ^ - A = 0 at z = 0. (4.1.20) 3z o

A normalized solution of these equations is:

cosh {K(h+z)}

o cosh (Kh) a (x,y), (4.1.21)

provided that the variable dimensionless wave number K is the real root of the dispersion relation

1 = K tanh (<h). (4.1 .22)

The two-dimensional amplitude function a (x,y) can be computed with Equation (4.1. 14), giving:

V . ( ^ VS ) = o (4.1 .23)

2 o K

O 2

/- n a because of the identity / A ^ dz =

J o 2 -h

I r , 2<h 1 with n = - {1 + . , .. , .J .

2 s m h (2<h)

Summarizing,the two refraction equations, written in dimensional form are: The "eiconal" equation defining the phase function S (x,y)

(^)^.(^)^ = V (^-'-^^^

with k the variable real wave number, depending on the variable local water depth h by the dispersion relation

m^ = gk tanh (k h ) . (4.1.25) o o

The two-dimensional energy flux equation, defining the amplitude function a (x,y),

„ n a ^ as . n a ^ as

#- (

°- ^ )

+ f

(-^

^ ) = 0, (4.1.26)

3x j^ 2 ax 3y ,^ 2

ay

o o in which 2k h n = i {1 + • K^i, t,J (4.1.27) 2 sinh (2k h) o

is a variable shoaling coefficient.

The eiconal Equation (4.1.24) is a well-known equation in the theory of geome-trical optics, defining the path of a wave ray.

Equation (4.1.26) is the energy flux equation from which can be derived the law of conservation of energy flux between two wave rays. This can be shown by apply-ing the divergence theorem of Gauss to the closed contour ABB'A'A of Figure 4.1.

wave orthogonals (rays)

wave fronts

s=constant

Figure 4.1 Definition Sketch Wave Rays and Wave Fronts

According to the expression (3.2.5) for the integrated mean energy flux vector P , this vector can be written as:

P * = l p ^ n a ^ VS

(4.1.28)

Along a wave ray, defined as the orthogonal of the line of equal phase S (see Figure 4.1), there must hold:

(m . VS ) = 0 ,

so the divergence theorem of Gauss gives the result, using Equation (4.1.26):

(4.1.29)

ƒ (m . P *) dn = ƒ (m' . P *) dn'

1 1'

which means that the energy flux between the two wave rays is constant.

The partial differential equation (4.1.24) can be transformed into a set of ordi-nary first order differential equations with the aid of the characteristic

direc-tion of the equadirec-tion. This characteristic direcdirec-tion appears to be the direcdirec-tion of a wave ray, which is defined as the orthogonal of a line of equal phase (see Figure 4.1).

The normal direction n at a line of equal phase (S = constant) is given by:

I

,3S 3S,

' '•Sx '

3y''-/

'-3X-'

^Sy'

Omitting the subscript zero of the asymptotic approximation for the sake of sim-plicity in notation, and putting (p,q) = (TT— , -TT—) and using the eiconal equation

ox dy

E " ^ (p,q). o

The coordinates of the orthogonal (wave ray) as a function of the distance E, along the ray are now given by:

and (4.1.30) (4.1.31) dx _ p dC " k dy ^ _q_ d?

k

Also there must hold

dp _ 9j) dx 9p dy

dC ~ 3x df" "^ ïy dë

dp 1 , 8p ^ 3p,

o •'

^, , ^. 3p 3q , 3'S 3^S , or, using the relation K ^ = -5^ (^-3— = ^ •^ ) dy dx dydx dxdy

dp 1 , 3p , 3q,

df = IT (P 3^ " ^ af)

o

which can be written with the aid of the eiconal equation as 3k

^ = ^ (4 1 32)

In the same way there can be derived: 3k

^ = — ^ . (4 1 33)

dC 3y (.t.i.Jj; The differential Equations (4.1.30)-(4.1.33) define the path of the wave ray.