Hamiltonian Formulation of Water Waves: 1D-formulation, numerical evaluations and examples

Hamiltonian formulation of

water waves

I D - f o r m u l a t i o i i , n t u n e r i c a l evaluations a n d examples

A . K . O t t a and M . W . Dingemans

Hamiltonian fornnulation H782 May 1994

Executive's summary

T h i s r e p o r t describes the f o r m u l a t i o n , n u m e r i c a l i m p l e m e n t a t i o n a n d a p p l i c a t i o n of a w e a k l y nonlinear wave m o d e l f o r finite d e p t h based o n a H a m i l t o n i a n f o r m u l a t i o n (see Radder, 1992). Due t o the t y p e of n o n l i n e a r i t y e x p l i c i t l y accounted f o r i n t h e expansion of the k e r n e l o f t h e H a m i l t o n i a n density ( s u m o f k i n e t i c a n d p o t e n t i a l energy per u n i t surface area), the m o d e l is v a l i d f o r waves of s m a l l , b u t finite a m p l i t u d e a n d f a i r l y l o n g wave l e n g t h (compared t o the water d e p t h ) i n r o u g h l y the same sense m a n y 'Boussinesq-type' models are. There are, however, a f e w significant differences. F i r s t l y , the H a m i l t o n i a n density i n the present f o r m u l a t i o n is always positive definite: a c o n d i t i o n necessary t o ensure good d y n a m i c a l b e h a v i o u r of the m o d e l equations f o r n u m e r i c a l c o m p u t a t i o n . Secondly, the dispersion e q u a t i o n o b t a i n e d f r o m t h e linearised version of the equations is exact. T h i s p r o p e r t y results i n a b e t t e r m o d e l l i n g of t h e phase relations (hence the wave a s y m m e t r y at a given l o c a t i o n ) of the s u p e r h a r m o n i c field w h i c h evolves f r o m the p r i m a r y wave system due t o surface n o n l i n e a r i t y . M o r e i m p o r t a n t l y , i t is possible t o remove the r e s t r i c t i o n o f waves being l o n g t h r o u g h a proper i n c l u s i o n of t h e 'short-wave' n o n l i n e a r i t y i n the expansion of t h e k e r n e l . T h i s results i n a u n i f o r m l y v a l i d m o d e l u n l i k e m o s t of the w e a k l y nonlinear models w h i c h are v a l i d over either deep or shallow water. F u r t h e r , i t is discussed i n the t e x t t h a t even f o r l o n g waves t h e 'short-wave' n o n l i n e a r i t y becomes l o c a l l y i m p o r t a n t near t h e crest of a wave as t h e surface curvature increases. I m p l e m e n t a t i o n of 'short-wave' n o n l i n e a r i t y is therefore considered as one of the first p r i o r i t i e s i n the f u t u r e developments o f t h e m o d e l .

T w o n u m e r i c a l models have been developed: a t i m e - d o m a i n m o d e l a n d a pseudo-spectral m o d e l based o n the sinc-series f o r the g l o b a l a p p r o x i m a t i o n s . T h e n u m e r i c a l code based on t h e sinc-series requires less c o m p u t i n g t i m e and gives the o p t i o n of choosing higher-order i n t e r p o l a t i o n f o r c o m p u t i n g derivatives and integrals over l o c a l i n t e r v a l s . We present here a figure as an example of t h e m o d e l p r e d i c t i o n of nonlinear e v o l u t i o n of a t r a i n o f non-breaking waves passing over an underwater bar. A n i n c i d e n t t r a i n of sinusoidal waves has been represented i n the c o m p u t a t i o n b y a packet of sinusoidal waves of finite l e n g t h w i t h the l e a d i n g edge a f e w wave-lengths b e h i n d t h e bar. D e t a i l s of t h e geometry a n d significance of t h i s test can be f o u n d i n chapter 7. I n

21 m

spite of a p r a c t i c a l disadvantage due t o t h e w a y the i n p u t conditions m u s t be specified i n t h e code (we have an i n i t i a l - v a l u e p r o b l e m ) these models can n o w be used t o s t u d y nonlinear evolutions of n o n b r e a k i n g waves over v a r y i n g d e p t h .

A t t e m p t s t o i n t r o d u c e effects of wave b r e a k i n g have n o t m e t w i t h m u c h success yet. I n

the f i r s t phase, a n i n t e g r a l c r i t e r i o n has been i m p l e m e n t e d t o determine i f t h e instanta-neous surface shape should lead t o b r e a k i n g . A p p l i c a t i o n of the c r i t e r i o n t o c o m p u t e d surface elevation f o r conditions observed t o have given rise t o m i l d b r e a k i n g i n labora-t o r y labora-teslabora-ts shows labora-t h a labora-t b r e a k i n g slabora-tage is n o labora-t reached. I labora-t is believed labora-t h a labora-t labora-t h i s difference between the c o m p u t e d results and the l a b o r a t o r y observations is caused b y the omis-sion o f t h e 'short-wave' n o n l i n e a r i t y ( u n d e r p r e d i c t i o n o f t h e surface steepness), r a t h e r t h a n the f a i l u r e o f the i n t e g r a l c r i t e r i o n . T h i s is another aspect w h i c h underscores the i m p o r t a n c e o f i m p l e m e n t a t i o n of 'short-wave' nonlinearity.

F i n a l l y t o conclude, comparison of t h e computed results w i t h the e x p e r i m e n t a l mea-surements and i n a w i d e r context w i t h several 'Boussinesq-type' models (see Dinge-m a n s , 1994a) provides strong Dinge-m o t i v a t i o n f o r f u r t h e r developDinge-ments o f t h e Dinge-m o d e l . T h i s is f u r t h e r s u p p o r t e d b y t h e inherent t h e o r e t i c a l appeal of t h e f o r m u l a t i o n s used. As summarised i n chapter 9, the m a j o r recommendations f o r f u t u r e w o r k are:

• T h e f o r m u l a t i o n of a b o u n d a r y v a l u e p r o b l e m instead of an i n i t i a l v a l u e p r o b -l e m . T h i s great-ly enhances the p r a c t i c a -l usefu-lness o f t h e p r o g r a m .

e I n c l u s i o n of short-wave nonlinearity. For a p p l i c a b i l i t y over the f u l l range of deep t o shallow water of a nonlinear wave m o d e l , i n c l u s i o n of short-wave n o n l i n e a r i t y is needed. F u r t h e r m o r e , i t has been shown t h a t f o r wave b r e a k i n g also the i n c l u s i o n of short-wave n o n l i n e a r i t y is needed t o o b t a i n m o r e realistic breaker indices.

• F u r t h e r s t u d y of wave b r e a k i n g characteristics a f t e r i n c l u s i o n of short-wave n o n l i n e a r i t y .

Hamiltonian formulation H782 May 1994

Contents

E x e c u t i v e ' s s u m m a r y i i 1 I n t r o d u c t i o n 1 2 T h e m a t h e m a t i c a l m o d e l 2 2.1 T h e H a m i l t o n i a n f o r m u l a t i o n 2 2.2 I D - f o r m u l a t i o n 3 2.2.1 S o l u t i o n i n the l ¥ - p l a n e 4 2.2.2 T h e inverse t r a n s f o r m a t i o n 4 2.2.3 Non-linear i n t e g r a l e q u a t i o n 6 2.2.4 L i m i t i n g cases 8 2.3 T h e e v o l u t i o n equations 9

2.3.1 E q u a t i o n s i n the physical plane 9 2.3.2 E x a c t e v o l u t i o n equations 12

3 N u m e r i c a l e v a l u a t i o n o f t h e t i m e - d o m a i n m o d e l 1 5

3.1 Regularised e v o l u t i o n equations 15

3.2 T i m e i n t e g r a t i o n 17 3.3 T h e choice of step size 18

4 S p e c t r a l a p p r o a c h 2 0 4.1 M e t h o d A : a;-space 20 4.2 M e t h o d A : x-space 22 4.2.1 E v o l u t i o n r a t e at fixed x-grids 22 4.3 M e t h o d s 23 4.4 M e t h o d C 24 4.4.1 E v o l u t i o n equation: r] 24 4.4.2 E v o l u t i o n equation: (p 25 4.5 S u m m a r y 26 5 N u m e r i c a l m e t h o d u s i n g s i n c - s e r i e s 2 7 5.1 N u m e r i c a l procedure 27 5.1.1 I n t e g r a t i o n procedure 27 5.1.2 C o m p u t a t i o n of the derivative 29 5.1.3 E v o l u t i o n rates 29 5.1.4 E n d t r e a t m e n t 30 5.2 T h e choice of step size 31

6 T h e p r o g r a m s 3 2

6.1 T h e m a i n p r o g r a m s 32 6.2 T h e i n p u t - g e n e r a t i n g programs solphi and p e r p h i 33

6.3 Periodic wave i n p u t 34 6.3.1 I n t r o d u c t i o n 35 6.3.2 T i e taper f u n c t i o n ƒ ( » ) 35 6.3.3 E x a m p l e of i n p u t file 37 7 N u m e r i c a l e x a m p l e s 3 9 7.1 S o l i t a r y waves 39 7.2 Periodic waves 45 7.3 C o n f r o n t a t i o n w i t i measurements 49 8 M o d e l l i n g o f b r e a k i n g 5 1 8.1 C r i t e r i o n f o r b r e a k i n g 51 8.1.1 B e i a v i o u r o f t i e i n t e g r a n d near g = 0 52 8.1.2 I n t e g r a l c r i t e r i o n i n p i y s i c a l space 52 8.2 M o d i f i c a t i o n of surface e v o l u t i o n 53 8.3 M o d i f i c a t i o n of v e l o c i t y e v o l u t i o n 54 8.4 S t u d y of b r e a k i n g c r i t e r i o n 55 8.4.1 V e r i f i c a t i o n of c o m p u t a t i o n of b r e a k i n g i n d e x 55

8.4.2 B r e a k i n g i n d e x f o r waves over a bar 56

8.5 S u m m a r y 57

9 C o n c l u s i o n s 6 1

Hamiltonian formulation H782 IVIay 1994

1 Introduction

I n t h i s r e p o r t t h e s t u d y i n t o the properties a n d n u m e r i c a l behaviour of t h e equations proposed b y Radder (1992) is described.

A t first, i n C h a p t e r 2, the m a t h e m a t i c a l f o r m u l a t i o n is given. T h i s follows closely Radder (1992). I t is i m p o r t a n t t o note t h a t t h e present t r e a t m e n t is s t r i c t l y v a l i d o n l y f o r I D wave p r o p a g a t i o n . T h i s is due t o t h e f a c t t h a t the methods rest u p o n c o n f o r m a l t r a n s f o r m a t i o n s . F u r t h e r m o r e , i n course of t h e d e r i v a t i o n , use is made of t h e a s s u m p t i o n t h a t the v a r i a t i o n of t h e s t i l l - w a t e r d e p t h is absent i n t h e H a m i l t o n i a n . T h i s m a y be compared w i t h the way i n w h i c h the mild-slope e q u a t i o n is derived: f o r t h e classical m i l d - s l o p e e q u a t i o n no v a r i a t i o n of t h e d e p t h is i n c l u d e d i n the H a m i l t o n i a n , b u t f o r t h e extended e q u a t i o n variations of t h e d e p t h i n the H a m i l t o n i a n have t o be i n c l u d e d , see Dingemans (1994b, Chapter 3 ) .

T h e n u m e r i c a l e v a l u a t i o n of t h e e v o l u t i o n equations f o r the t i m e - d o m a i n m o d e l is described i n chapter 3.

A spectral approach based o n an expansion i n terms of sine f u n c t i o n s is described i n chapter 4. T h e basic idea b e h i n d such a spectral description is t o o b t a i n i m p r o v e d n u m e r i c a l elRciency. T h e t i m e - d o m a i n m o d e l s t i l l takes m u c h computer t i m e , especially because t h e e v a l u a t i o n o f t h e kernel l o g t a n h ( ) w h i c h has t o be evaluated over t h e w h o l e c o m p u t a t i o n a l d o m a i n . E a r l i e r a t t e m p t s t o t e r m i n a t e t h e c o m p u t a t i o n o f t h i s k e r n e l b e y o n d a l o c a l i n t e r v a l t o take advantage o f the r a p i d decay of t h e k e r n e l h a d no success. I n s t a b i l i t i e s developed t h e n v e r y soon f r o m near the ends of t h e geometry. I t should be emphasised t h a t spectral s o l u t i o n based o n t h e Fourier-series expansion over a finite span has been f o u n d t o be unsuccessful f o r t h e present H a m i l t o n i a n f o r m u l a t i o n . T h r e e possible approaches based o n sine series have been discussed i n chapter 4. T h e n u m e r i c a l procedure f o r one of the three approaches has been described i n chapter 5.

A short d e s c r i p t i o n of t h e programs is given i n chapter 6. Some n u m e r i c a l examples have been given i n chapter 7.

A n approach t o m o d e l l i n g wave b r e a k i n g has been discussed i n chapter 8, f o l l o w e d b y the conclusions i n chapter 9.

T h e research has been carried o u t b y A . K . O t t a a n d M . W . D i n g e m a n s , w h i l e t h e constant guidance o f A . C . Radder of R i j k s w a t e r s t a a t has been i n v a l u a b l e .

The mathematical model

2.1 T h e Hamiltonian formulation

T h e governing equations f o r i r r o t a t i o n a l wave m o t i o n o n an incompressible i n v i s c i d fluid are given b y t h e Laplace equation and the three b o u n d a r y conditions:

= 0 -h{x) <z< Cix,t) 0 $ dt — + V $ - V / i = 0 dz z^Ci^,t) z = C{x,t) z = —h(x), (2.1a) (2.1b) (2.1c) (2.1d) where

T

[C]

=

7 V VC

(l + IVCP)

1/2

(2.1e)

a n d 7 is the surface tension a n d two-dimensional vectors have been used, x = {x,y^. I t has been s h o w n t h a t the H a m i l t o n i a n constitutes a v a r i a t i o n a l p r i n c i p l e w h e n i t is expressed i n t e r m s of t h e free surface elevation C a n d the value o f t h e velocity p o t e n t i a l at t h e free surface ( f { x , t ) = ^{(x,({x,t),t}, see, e.g., Zakharov (1968) a n d Broer (1974). T h e t o t a l energy of t h e fluid is given b y

n = J J dxdyH = J J dxdy (V + T )

< J J dxdy I^ K ^ + 7 ^ 1 + |VCP - 1

+

_{/ / 4}

(9$ ( 2 . 2 )

I n the sequel of t h i s r e p o r t no effect of surface tension is accounted f o r , i.e., T = 0 is t a k e n . T h e n we have

V=^pgC' a n d T =\p J^_^ dz (2.3)

T h e v a r i a t i o n of t h e H a m i l t o n i a n , given b y (2.2), is equivalent w i t h the o r i g i n a l f u l l set of water wave equations. However, one has the o p p o r t u n i t y t o find a p p r o x i m a t e equat i o n s f r o m an a p p r o x i m a equat i o n of equathe H a m i l equat o n i a n w i equat h equathe advanequatage equat h a equat good d y n a m

-Hamiltonian formulation H782 May 1994

i c a l b e h a v i o u r is guaranteed w h e n t h e p o s i t i v e definiteness of t h e exact H a m i l t o n i a n density H is also c a r r i e d over t o t h e a p p r o x i m a t e density Ha- E v e r y a p p r o x i m a t i o n Ha w h i c h m a i n t a i n s t h e symmetries of H guarantees a u t o m a t i c a l l y t h e corresponding conservation laws t o h o l d .

N o t i c e t h a t t h e H a m i l t o n i a n density is effectively a f u n c t i o n of t h e free surface variables ( a n d (p. V a r i a t i o n o f t o { a n d ip gives t h e t w o canonical equations:

SH

dc .

SH

dip

where S denotes the v a r i a t i o n a l derivative. N o t i c e t h a t Eqs. (2.4) are t h e free sur-face conditions ( 2 . 1 b ) a n d (2.1c) w h i l e Eqs. (2.1a) and ( 2 . I d ) have been used as side

c o n d i t i o n s . T h e m a j o r d i f f i c u l t y i n a p p l y i n g t h e H a m i l t o n i a n is t h a t t h e k i n e t i c en-ergy density should be f o r m u l a t e d i n t e r m s of t h e surface elevation ( a n d t h e v e l o c i t y p o t e n t i a l (p along t h e free surface.

I n t h e sequel we w i l l use a H a m i l t o n i a n w h i c h is scaled w i t h t h e density: H := Hip a n d s i m i l a r l y scaled p o t e n t i a l a n d k i n e t i c energy ( V a n d T respectively) t o c o m p l y w i t h the n o t a t i o n used i n Radder (1992).

2.2 ID-formulation

For problems w i t h one h o r i z o n t a l dimension an a p p r o x i m a t e H a m i l t o n i a n can be de-veloped w h i c h is u n i f o r m l y v a l i d f r o m deep t o shallow water ( R a d d e r , 1992). I n t h i s case, the k i n e t i c energy d e n s i t y T can be w r i t t e n as

2^ dx 2dx^'

where i/) = ^yp^ {x,({x,t),t} is t h e s t r e a m f u n c t i o n evaluated at t h e free surface. (2.5)

To express i/> i n t e r m s of t h e canonical variables (p a n d a c o n f o r m a l m a p p i n g o f t h e f l u i d d o m a i n Z = x + iz i n t o an i n f i n i t e s t r i p i n t h e complex T1^—plane W = x + "^C '^^ used (Woods, 1961)

C(x') + c o t h TT

{w-x')

Kx') .

(2.6)

A s o l u t i o n f o r t h e s t r e a m f u n c t i o n ^{x,z,t) at t h e free surface is sought f o r i n t w o steps:

1. Solve t h e p r o b l e m i n t h e 11^—plane

2. F i n d t h e inverse t r a n s f o r m a t i o n %(a;) along t h e free surface.

2.2.1 S o l u t i o n in t h e W-p\ane

A s o l u t i o n of t h e linear p r o b l e m ^ + = O i n O < ^ < l a n d the b o u n d a r y conditions ^ = 0 at = 0 and $ = '4>{x) at ^ = 1 b y means of Fourier t r a n s f o r m s gives at the free surface ^ = 1

where l o g is the n a t u r a l l o g a r i t h m . I n the t r a n s f o r m e d plane t h e H a m i l t o n i a n H becomes

a n d w i t h V>(x) s u b s t i t u t e d f r o m (2.7) one obtains

w i t h t h e kernel TZ{x-,x') given b y

7e(x,x')

= - - l o g t a n h

Ix - x'|) •

(2-10)

2.2.2 T h e inverse t r a n s f o r m a t i o n

To express (2.7) i n terms of variables of the p h y s i c a l plane, an expression f o r the f u n c t i o n x(a^) along the free surface has t o be f o u n d . T h e i m a g i n a r y p a r t o f the t r a n s f o r m a t i o n (2.6) reads

- ^ / I c o s l . ( x - ° x 1 - c o . . { " ' ^ ' " ' ^ ' <^-"'

T h r o u g h a Fourier t r a n s f o r m of t h i s i m a g i n a r y p a r t a symbolic operator e q u a t i o n is o b t a i n e d ( U a d d e r , 1992, E q . ( 2 6 ) ) :

/ 1 ^ ( I \

T h i s e q u a t i o n has t o be i n v e r t e d t o o b t a i n t h e required f u n c t i o n

xi.^)-

T h e subsequent analysis presented b y Radder is v a l i d i n the s t r i c t sense f o r t h e s i m p l i f i e d case o f a h o r i z o n t a l b o t t o m f o r w h i c h (2.12) simplifies t o

Hamiltonian formulation H732 May 1994

dx__d_ r_d_

dx dx \ \dx

vix)

w h e r e t h e n o t a t i o n r]ix) =

h + C(x)

(2.13a) (2.13b) is used. A n unlcnown f u n c t i o n s(x) is i n t r o d u c e d t h r o u g h dx

; ^ = ( l

+ e ) , . ( 2 . 1 4 )

N o t i c e t h a t t h i s r e l a t i o n at t h i s stage is m e r e l y another w a y of w r i t i n g r e l a t i o n (2.13a).

T h e k i n e t i c energy T can be r e w r i t t e n i n t e r m s o f t h e physical variables i n t h e f o l l o w i n g way. F i r s t notice t h a t d dx'^^ and dip^dipdx^ dx

dx

dx

dx '

so t h a t ld(f f , . d ( p ^ , 1 dx

f,,dx'

d x ' , z^T =^ J dx(px J dx'(p^iR{x,x')

I n t r o d u c i n g the k i n e t i c energy density b y

dx'(pa,(px''^eix,x') , (2.15)

we t h e n have T = f dxT . T h e kernel TZ^ now follows f r o m (2.10) b y n o t i n g t h a t (2.14) gives

dx =

dx W e t h e n have t h e s y m m e t r i c f u n c t i o n TZ^ given b y : / / ^ 11 ^ T, /^'^ r ' IZs{x,x;7]) = l o g t a n h — / TT V 4 Ja, dr

x - x =

_I

_{( 1} dr (2.16)

I t is n o t e d t h a t t j i e k i n e t i c energy f u n c t i o n a l , expressed i n the x-variable, is positive definite, a n d t h i s p r o p e r t y is preserved o n t r a n s f o r m i n g f r o m x t o x w h e n the J a c o b i a n (2.13a) is positive and b o u n d e d .

T h e H a m i l t o n i a n i n the physical plane t h e n is

^^ = 7, dxigCi- dx'cpa:<pa;'ne{x,x';ri)) . (2.17)

^ J — CO \ J — (X> J

E i t h e r t o be able t o use (2.17) d i r e c t l y t o o b t a i n e v o l u t i o n equations f r o m , or t o trans-f o r m results obtained i n the % plane t o the physical plane, i t is necessary t o trans-f i n d a p r e s c r i p t i o n f o r t h e u n k n o w n f u n c t i o n e{x). T h i s is t h e subject of next subsection.

2 . 2 . 3 N o n - l i n e a r i n t e g r a l e q u a t i o n

We n o w seek an e x p l i c i t expression f o r t h e f u n c t i o n £ ( x ) . W e f i r s t notice t h a t equation (2.13a) can be w r i t t e n i n the f o r m of a nonlinear i n t e g r a l e q u a t i o n (Radder, 1992) as

^ ^ - ï ï U J ^ ^ ] _ ^ sinh^iTTx' ^ (2.18)

To solve t h i s nonlinear i n t e g r a l equation, an expansion i n p a r t i a l f r a c t i o n s , w h i c h is v a l i d f o r a l l values of the argument, is used b y Radder (1992):

dx

OO I m-1 m-1 + — mTT 4-. + m-K dx V • (2.19) W i t h the n o t a t i o n rx dr P

Jo

V ^

dp

^

dx

' (2.20)

and the f u n c t i o n e{x) defined according t o (2.14), an e q u a t i o n f o r e is obtained f r o m the expansion (2.19) as

m = l

1

+

V (2.21)

Xl + s)D-nnr (1 + s) D + mir

To solve t h i s equation f o r e{x), Radder introduces the operators G^™^ and GQ™^ b y

D , „r™) D = and GS™) = ^

and because Gi"^^ can be expressed i t e r a t i v e l y i n terms of GQ™^ we have

(2.22)

) ^ G ( ™ ) ( £ G ( ' " ) ) ' .

For t h e f u n c t i o n e{x) we t h e n o b t a i n the equation

G'(-) = ^ ( _ l ) A G M ( £ G M y

A=0

Hamiltonian formulation H782 May 1994

A=0 w i t h

-| OO r 1,

V

E q u a t i o n (2.21) or (2.24a) can now be solved i t e r a t i v e l y as

Io Io - I l £0 = 0 , £ l =

£2

l - I o ' l - I o + I i

where Iq is given t o be (see Radder, 1992)

dq d Io

V(p) Jo

Tjip) Jo exp(7rg) -Idq a n d II follows as Hp +

l) + v{p

-1)] , (2.24a) (2.24b) (2.25) (2.26) 1 _

j"^ r

[p] Jo Jo

dqdq"

vip)

dq exp TT (g + g') - 1

£1

(p +

e)-^vip

+ q + q') + £i{p- q)-^v{p-q-a')

(2.27) S o l u t i o n f o r £fc, ^ = 0 , 1 , 2 , • • • can i n p r i n c i p l e be o b t a i n e d t h r o u g h a r e c u r r i n g proce-dure i n t h i s way.

F i n i t e d e p t h

For water of f i n i t e d e p t h we have-"^ \\D\\ < TT a n d a T a y l o r series expansion of GQ™^ m a y be used (Radder, 1992, E q . ( 4 3 ) ) :

mir

rmr

\m7rj

yrmr J

(2.28)

a n d liliewise f o r "^K T h e expansion of (2.24a) m a y t h e n be given i n a more conve-n i e conve-n t f o r m . T h e f o l l o w i conve-n g quaconve-ntities are i conve-n t r o d u c e d :

Fk = -D''7j; Fo = l . V

T h e f o l l o w i n g recurrence relations m a y t h e n be derived f r o m (2.29) a n d (2.20):

F (r^Fk) = i]Fk+i a n d DF^ = F^+r - Ff,Fi .

(2.29)

(2.30)

^A bounded opeiatoi T is defined by (ƒ, V) < (ƒ, ƒ) with f some well-behaved function and (•, •) an innei pioduct. The norm of the opeiatoi, HT'H is defined as the maximum of ( V f , V f ) /(ƒ, ƒ).

SubstitTition of t l i e series expansion (2.28) i n (2.23) and using the recurrence relations (2.30) t h e n leads t o an expansion of e i n terms of the f u n c t i o n s F},. U p t o f o u r t h order Radder obtains (his e q u a t i o n ( 4 6 ) ) :

£ = -^i^2 - ^ (I4 - - lOFi + bF2Fl) , (2.31)

w h i c h can be r e w r i t t e n as^ (cf. Radder, 1992, E q . ( B 6 ) )

{-^vt - 6 (4r;^ + VV<v) mxx + '^''if'Vxx'nxxx + V^Vio) — •

(2.32)

I t is i m p o r t a n t t o note t h a t crucial f o r t h i s a p p r o x i m a t i o n is the v a l i d i t y of the ex-pansion ( 2 . 2 8 ) . I n t h e f i r s t place t h e exex-pansion is o n l y v a l i d f o r water of f i n i t e - d e p t h , and i n the second place the expansion is a T a y l o r t y p e expansion, n o t y i e l d i n g p o s i t i v e -d e f i n i t e a p p r o x i m a t i o n s f o r the class of square integrable f u n c t i o n s w i t h ?/ > 0. Perhaps a better a p p r o x i m a t i o n can be o b t a i n e d b y using a Fade - t y p e a p p r o x i m a t i o n f o r t h e f u n c t i o n s g'q^\

W i t h the H a m i l t o n i a n given b y (2.17), the e v o l u t i o n equations f o r the canonical v a r i -ables C and are i n p r i n c i p l e given b y the defining equations (2.4). However, due t o t h e c o m p l i c a t e d dependence o n £, d e r i v a t i o n o f the e v o l u t i o n equations is n o t a t r i v i a l task, except f o r the case o f £ = 0.

2.2.4 L i m i t i n g cases

Consider the expression (2.16) f o r TZ^. For a constant d e p t h o f h, we can w r i t e

7le{x, x'; rj) = l o g t a n h _{( - - r '} V 4

hJx

dr

(! + £)(! +CM)

(2.33)

L i n e a r t h e o r y is o b t a i n e d f o r £ = 0 a n d 77 = h (i.e., i n f i n i t e s i m a l waves, a n d t h u s ( is neglected); the expression f o r Re = Ri t h e n coincides w i t h the one given b y B r o e r (1974) f o r linear waves. For deep water i t t u r n s out t h a t s = 0{ka) a n d t h e n Stokes' t h e o r y is a p p r o p r i a t e , rj = h, so t h a t

Re{x,x';h) = l o g t a n h Ah

dr

+ £

(2.34)

For shallow water we have e = 0{kh • ka) and we can p u t £ = 0 (Boussinesq a p p r o x i -m a t i o n ) . T h e n :

^With Fk = {v~^D''ri) /v = [iv^xf] V we have Fi = {l/r])Tida:Ti = and F2 = [(1/??) {vda:f] V = dm{r)dx)r) and thus, F2 = vl + VVxx- For Fs we have ƒ3 = {Tjdx (ndx)) V, yielding F3 =

Hamiltonian formulation H782 May 1994 Ro{x,x';r)) = - i l o g t a n h ( ^

L

^' dr 1 + ^ (2.35)

T h e advantage o f a H a m i l t o n i a n m o d e l based o n (2.35) over a Boussinesq-type m o d e l is t h a t i n t h e f o r m e r m o d e l the linear t h e o r y is recovered i n case of deep w a t e r . For p r a c t i c a l applications, such as t h e e v o l u t i o n of a s o l i t a r y wave over an uneven b o t t o m , n u m e r i c a l s o l u t i o n o f t h e r e s u l t i n g i n t e g r a l equations is needed (see Z w a r t k r u i s , 1991).

2.3 T h e evolution equations

I n t h e f o l l o w i n g we describe t h e e v o l u t i o n equations i n t h e p h y s i c a l plane i n c l u d i n g the short-wave n o n l i n e a r i t y f u n c t i o n e w h i c h depends o n rj a n d i t s derivatives. T h e r e a f t e r , we l o o k at the possibihty of the e v o l u t i o n equations i n the t r a n s f o r m e d plane.

2.3.1 E q u a t i o n s in t h e physical plane

S u r f a c e e l e v a t i o n C

To o b t a i n t h e e v o l u t i o n equations i t is necessary t o evaluate the v a r i a t i o n a l derivatives. U s i n g the general f o r m u l a t i o n

sn

dH d fdH 6ip d(p dx \dipx we n o t e f o r t h e H a m i l t o n i a n H, given b y ( 2 . 1 7 ) , t h a t (2.36) dt

' èif

d (dE d f dH dx \dipxJ dx' \d(p^ = ƒ dx''^Re{x,x'-rj{x)) dx^Re{x,x';v{x)) = - J dx'cpxi^ - J dx<px^ . (2.37) D u e t o t h e s y m m e t r y of Rs{x,x';ri) w i t h respect t o x and x', (2.37) can also be

ex-pressed as dC dt = ƒ dx'ipa;'Re{x,x']r}{x)) = - J dx'ip.. dRe ^' dx (2.38) T h e d i f f e r e n t i a t i o n o f R^ proceeds as f o l l o w s . W e i n t r o d u c e t h e n o t a t i o n 1' dr

+

delft hydraulics 9

1 TT F r o m Re given i n (2.16) we have 71^ = l o g t a n h —v a n d t h u s , TT 4 dRe _ 1 dvldx dx 2 sinh ( f t ; ) " T x r - ^ dv . . ^.dê dê 1 W i t h — = sign('!7)-— a n d — = . , , we o b t a i n dx ^ ^ Ux dx ( l + £ ) ? ; ( a ; ) dRe _ 1 sign(i?) dx 2 ( l + £)7?(a;) s i n h ( f - y ) " Consider n o w t h e t e r m sign('!?)/sinh ( f t ; ) : sign(i?) sign('!?) + 1 (2.39)

sinh ( f - y ) sinh ( f | i ? | ) sinh ( f i ? )

- 1 s i n h ( - f i ? ) s i n h ( f i ? ) i f 1? > 0 i f 1? < 0 . We can therefore w r i t e dRe 1 1

dx 2 ( 1 + £) ri{x) sinh ( f i?)

1 1

(2.40) n i + e),ix) s i n h d r ^ ^ )

T h e e v o l u t i o n e q u a t i o n f o r t h e free surface elevation ( f i n a l l y becomes

V e l o c i t y p o t e n t i a l </?

T h e e v o l u t i o n e q u a t i o n f o r the free surface wave p o t e n t i a l ip is m o r e d i f l i c u l t t o o b t a i n . T h i s is due t o t h e f a c t t h a t e{x,t) is a f u n c t i o n of rj{x,t) a n d therefore also a f u n c t i o n of

C(*?*)-

Since h is o n l y a variable constant,

6r] = Sh + 6C = SC . (2.42) T h u s , 6( m a y be replaced b y 6r]. N e x t , SH, the increment i n t h e H a m i l t o n i a n due t o

^77, a s m a l l i n c r e m e n t i n t] and the associated Ss can be expressed as

Hamiltonian formulation H782 May 1994

/

CO -OO

i

r

r

—

4

7 - c

Since

1^1

4 j _ o o (l+^)'7

8rj

St] Jx' 1 + {l + e)r] (2.43) — one lias f r o m (2.43) s i n l i f | Ö | s i n h f ö

/

OO gCSCdx' --OO 1 fOO POO - I dx' dx" / » ^ N

/

OO gCèCdx' + -OO

(r

'^(iT

OO s i n l i I ^x"Px"

ii:"

(1+^)77 dr (1 + e)

V

(2.44)

One m a y proceed f r o m here t o derive a suitable expression b y c a r r y i n g out i n t e g r a t i o n b y p a r t s o f the r i g h t - h a n d side. A l t e r n a t i v e l y , a similar approach as used i n appendix E of Radder (1992) m a y be adopted. Since (2.44) must h o l d f o r a l l 6r) G C , St] m a y be replaced b y any sequence fn(x) approaching t h e delta f u n c t i o n ^(a;).

T h e f o l l o w i n g properties r e l a t e d t o t h e delta f u n c t i o n w i l l be u s e f u l d u r i n g the deriva-t i o n . F i r s deriva-t , as i n Radder (1992, A p p e n d i x E ) , use w i l l be m a d e of

fx"

/ drF(r)S{r - x) = F{x)U{x - x')U{x" - x) Jx'

(2.45)

where F{x) represents a well-behaved f u n c t i o n and U{x) t h e u n i t Heaviside f u n c t i o n . W e f u r t h e r use t h e f a c t t h a t {e.g., Greenberg, 1978; p . 71)

rx"

/ drF{r)Sj{r - x) = {-iyFj{x)U{x - x')U{x" - x). Jx'

(2.46)

T h e subscript j i n (2.46) denotes t h e order of d i f f e r e n t i a t i o n of the f u n c t i o n . We n o w t u r n t o evaluate Se. U s i n g (2.32) f o r e and c a r r y i n g out the o p e r a t i o n t o 0(kay, we have Se(r) = + V^Vrr + VrrSri) • A p p r o x i m a t i n g Srj^r) b y S{r — x) i n (2.47), one has (2.47) Se = --[2r]rSi{r - x) + 7]S2{r - x) + VrrK'^ - x)] (2.48) Hence, delft liydraulics 11

Jx' (1

+ £)V

1"

Jx'

dr êrj

₊

Se

= U{x-x')U{x" -x)Bi{x) where B\ up t o 0 [(fca)^] is giveu b y

(2.49) 5 i ( 7 ? , £ ) ( l + £)7?2

+

d / 27;r7,

+

9 2 . ( 1 + 6 ) ^ 2 d x \ { i ^ e f r j ^ ) dx-^y^l^efri^,

b y using a reduced f o r m of e(a;); i.e.,

U s i n g (2.49) i n (2.44) yields (2.50) (2.51) SH

= K +

/

OO poo dx' / t^a;" •00 J—00 4 j-00 <Px"Px"

00 sinh I

0

a;" rfr Ï7(a; - x')U(x" - x) . (2.52)

R e w r i t i n g the i n t e g r a l o n the second l i n e of (2.52), we o b t a i n

Consequently, t h e r e s u l t i n g e v o l u t i o n equation f o r (p is

(2.53)

sn

1 fx fOO

- 9 C - M v , e ) dx' dx" ^Px^Px" (2.54)

I t m a y be n o t e d t h a t equation (72b) of Radder (1992) results as a special case of (2.54) f o r £ = 0.

2.3.2 Exact e v o l u t i o n e q u a t i o n s

A s discussed so f a r a p p r o x i m a t i o n s have t o be used f o r t h e inverse t r a n s f o r m a t i o n i n order t o express the e v o l u t i o n equations i n the physical plane. A n a l t e r n a t i v e is t o investigate t h e w o r k a b i l i t y of using the exact equation o n the t r a n s f o r m e d plane. T h e

Hamiltonian formulation H782 May 1994

exact H a m i l t o n i a n H a n d 7^(x,x') are given respectively b y (2.9) a n d (2.10).

S u r f a c e e l e v a t i o n C

For t h e e v a l u a t i o n o f SH/Scp we need o n l y t h e k i n e t i c energy p a r t a n d we have

6n_ST_

Since t h e variable x depends o n 7], we need t o r e w r i t e t h e exact expression of T p a r t l y i n terms o f t h e p h y s i c a l variable x t o be able t o o b t a i n t h e v a r i a t i o n a l derivative 8T /Sip.

Observing t h a t dxl^Px' = dx'fx'i we have

r=\jdxcpx jdx'<px'R{x\x;V) • (2-55)

. , 8r d f dT\ d f d r \

Since T has no e x p l i c i t dependence cp i t s e l f , we have — = " " ^ ^ \d(p ) ~ ~dx \dïp' a n d we o b t a i n

-f = - ^ dx'cp,>n{x,x') • (2-56)

8(p ox J-oo

A l o n g t h e free surface ^ = 1 we have -w^ = -t^tt- a n d therefore 8H/8ip can be e x a c t l y OX CIX

expressed as

dx (dx\~^

Since one has — = I —- subject t o t h e c o n d i t i o n t h a t t h e Jacobian is n o t zero dx \dxJ

( w h i c h is t r u e f o r t h e m a p p i n g used), dx/dx can be calculated exactly f r o m

(2.12)

or f o r a h o r i z o n t a l b o t t o m f r o m

(2.18).

T h u s ,

(2.57)

defines t h e e v o l u t i o n e q u a t i o n f o r ( e x a c t l y i n t h e t r a n s f o r m e d variables.

V e l o c i t y p o t e n t i a l <p

T h e exact e v o l u t i o n equation f o r (p f o l l o w s f r o m t h e free surface conditions t o be

where again a l l t h e derivatives can be expressed i n t h e t r a n s f o r m e d plane b y using t h e Jacobian o f t h e t r a n s f o r m a t i o n .

A n i m p o r t a n t p o i n t t o note is t l i a t the e v o l u t i o n rates, given b y (2.57) a n d (2.58), correspond t o t h e rates at a f i x e d physical l o c a t i o n x and n o t at a f i x e d %. A f i x e d X-node does n o t correspond t o the same x - l o c a t i o n due t o the change i n t h e surface elevation w i t h t i m e . For t h i s reason, the n u m e r i c a l values need t o be m a p p e d f r o m the set of x-nodes t o the set of the f i x e d x-locations and vice-versa at each time-step. T h i s is one o f t h e basic disadvantages of using the exact e v o l u t i o n equations o n the x-space.

Hamtlton'tan formulation H782 May 1994

3 Numerical evaluation o f t h e time-domain model

I n t l i e present chapter, we describe the n u m e r i c a l m o d e l w h i c h is based o n t h e f o r m u -l a t i o n s i n the p h y s i c a -l p-lane. T h e n u m e r i c a -l m o d e -l described here w i -l -l be -l a t e r referred t o as ' h a m i l T ' .

3.1 Regularised evolution equations

A s t h e integrands i n (2.38) and (2.54) are singular at a; = a;', regularised f o r m s are f i r s t derived i n a way s i m i l a r t o t h e one adopted i n Z w a r t k r u i s (1991). F u r t h e r , we derive a m o d i i i e d f o r m of (2.54) where t h e outer double integrals need n o t be n u m e r i c a l l y evaluated. S u r f a c e e l e v a t i o n ( W i t h dr ( l + £ ) r / we w r i t e

/

OO fOO d a ; ' A ( a ; ' , a ; ) l o g t a n h ! ƒ ( ? ] , £ ) | = / (ia;''y(a;')logtanh | / ( 7 7 , £ ) | --OO J — OO ^ T / X /•'^ , / l o g t a n h I f (7?,£) I , ^, K l + e ) , l w / _ ^ x ' ^ J | J i . (3,1) I t follows f r o m t h e d e f i n i t i o n of / ( T / , £) b y using L e i b n i t z ' s r u l e t h a t N o t i n g t h a t ƒ ^ oo as a;' —> oo a n d ƒ ^ — oo as — oo we can f u r t h e r w r i t e . o o ^ " ' [ ( l + £ ) ^ ] ( x O ^ ° ' * " ^ ' ^ ^ ' ' ' ' ^ ' = - y _ c ^ / l o g t a n h | / ( . , . 0 | = -TT . ( 3 . 3 ) delft hydraulics 15

U s i n g (3.3) i n (3.1) yields

/

OO da;'A(a;',cc)logtanh \f{x^x^) -OO T T K I + e)77](a;) +

/

OO da;''y(a;')logtanli|/(a;,a;')| . (3.4) -OO A l t e r n a t i v e l y ,

/:

dx'v{x')logta,niL\f(x,x')\ = - 7 r [ ' y ( l + £)r/](a;) roo + / rfa;'A(a;',a;)logtanh|/(a;,a;')| . (3.5)

T l i e i n t e g r a n d X{x', x) l o g t a n h | / ( s , a;')| is regular at a; = a;' a n d t l i e n u m e r i c a l integra-t i o n o f integra-t h e r i g h integra-t h a n d side o f (3.5) can be carried o u integra-t b y using a s integra-t a n d a r d n u m e r i c a l technique like t r a p e z o i d a l r u l e (or Gauss-quadrature f o r higher accuracy). T h u s , t h e s t a r t i n g e q u a t i o n f o r t h e n u m e r i c a l e v o l u t i o n o f t h e free surface is TT dx

/

oo rfa;'A(a;', a ; ) l o g t a n h I / ( a ; , « ' ) -OO (3.6) V e l o c i t y p o t e n t i a l (p

T h e double i n t e g r a l and t h e singular behaviour o f t h e denominator at x' — x" i n (2.54) make a direct n u m e r i c a l i m p l e m e n t a t i o n d i f f i c u l t . Expression (2.54) is f u r t h e r s i m p l i f i e d as f o l l o w s . DifFerentiation o f (2.54) w . r . t . t o x gives f i a gCx -1 „ d 1 ^ ] r dx' r dx"g{x',x") dx J J—oo Jx \ H /'OO --Bif-^ / dx' / dx"g{x',x") Z dX J—cx> Jx where g{x',x") f x ' f x " U s i n g Leibnitz's r u l e , one has

dx

/

X roo rx -7 /•oo dx' dx"g{x',x") = / dx'^ dx"g{x',x") •OO Jx J—OO dX J(jQ fOO + / dx"g{x,x") Jx (3.7) (3.8)

/

X roo dx'[-g{x',x)]+ / dx"g{x,x") ( 3 . 9 ) -OO J X

Hamiltonian formulation H782 May 1994 I t f o l l o w s f r o m t h e d e f i n i t i o n of g t h a t g{x',x) = -g{x,x'). Hence, (3.10) dx 7 nx rco rx poo - dx' dx"g{x',x") = / dx'[-g{x',x)]+ dx"g{x,x") X J—oo Jx J-oo Jx

/

OO dx'g{x,x'). (3.11) -OO W e r e c a l l f r o m (2.41) t h a t

/

°° dx'g(x,x') = -2{l+ e)wxCt. (3-12) -oo U s i n g (3.11) and (3.12) i n (3.7) we get f i x = -gCx + dBr j dx'{l + e)Wx'Ci J—oo dx + 5 i ( l + e)r]ipxCt. S u b s t i t u t i n g v f o r (px, (3.13) is f i n a l l y expressed as Vt = -gCx +

dx _J—oo r dx' [v{l + e)vCt] + Bl [v{l + e)vQ

(3.13)

(3.14)

Expression (3.14) is t h e s t a r t i n g e q u a t i o n f o r the n u m e r i c a l e v o l u t i o n of v (or conse-q u e n t l y o f (p). T h e regularised set o f econse-quations is t h e n given b y

Ct = - ¬

TT dx

Vt = -gCx +

7 r [ ' u ( H - e ) r / ] ( « ) + / (ia;'A(a;',a;)logtanh|/(a;,a;')

dx _J—oo r dx' [v{l + e)riCt] + Bl [v{l + e)v(t] •

(3.15a)

(3.15b)

I n t h e n u m e r i c a l code 'h a m i l T' , t h e derivatives i n the e v o l u t i o n equations are de-t e r m i n e d b y using a f o u r de-t h - o r d e r m i d - p o i n de-t rule^ a n d de-the i n de-t e g r a de-t i o n is carried oude-t according t o the t r a p e z o i d a l r u l e . Nodes are u n i f o r m l y d i s t r i b u t e d over t h e e n t i r e d o m a i n .

3.2 T i m e integration

T i m e i n t e g r a t i o n of t h e t w o e v o l u t i o n equations (3.15a) and (3.15b) are p e r f o r m e d according t o either an e x p l i c i t scheme ( E u l e r ) or a p r e d i c t o r c o r r e c t o r scheme ( A d a m s

-Lagrangian polynomial approximation based on discrete values at 5 nodes are used, with the node of interest being the middle one. Higher-order interpolation can also be later implemented.

B a s h f o r d - M o u l t o n , liencefortlx referred t o as ' A B M ' ) . F r o m k n o w n values of ( and v at a n i n s t a n t t, values at t h e next i n s t a n t At l a t e r are given i n t h e Euler procedure b y

C{x,t + At) = ((x,t) + A t ^ ^ ^ ^ + OiAty (3.16) v{x,t + At) = i ; ( x , i ) + A i ^ ^ ^ + 0 ( A t ) 2 (3.17)

where a n d Vt are c o m p u t e d f r o m t h e e v o l u t i o n equations (3.15a) a n d (3.15b).

For a n o r d i n a r y d i f f e r e n t i a l e q u a t i o n dy/dt = f ( y , t ) , t h e value of y At l a t e r is given b y a two-step procedure i n t h e A B M scheme:

= / + ^ (552/^° - 5 9 % - i + 37%-2 - 9%-3) ( g . i g a )

= ^ ° + i i^y''+^^y" - ^y^'+y^') (^-i^^)

where t h e superscripts denote t h e number of time-steps w i t h respect t o t h e present i n s t a n t , Ip b e i n g t h e p r e d i c t e d value a n d 1 t h e corrected one. I n t h e context of t h e present e v o l u t i o n equations, t h e p r e d i c t e d values of ( a n d v are evaluated f i r s t . T h e e v o l u t i o n rates (l^ a n d -yj^ are t h e n calculated f r o m (3.15a) o n t h e basis o f the cal-c u l a t e d C^ï' a n d v^P. Sincal-ce e v o l u t i o n rates f r o m three previous time-steps are used i n t h e A B M procedure, t h e e x p l i c i t Euler scheme is used f o r a f e w time-steps at t h e s t a r t of t h e c o m p u t a t i o n . For t h i s purpose, a m u c h smaller time-step ( f o u r times smaller) is used t h a n t h e i n t e n d e d time-step t o be eventually used w i t h t h e A B M scheme. W e n o t e t h a t t h e t r u n c a t i o n error is of O ( A t ^ ) i n t h e E u l e r scheme w h i l e t h a t of t h e A B M scheme is 0 (At^).

3.3 T h e choice of step size

T h e guidelines t o choose step sizes are as follows. T h e s p a t i a l step Aa; is chosen f i r s t based on some i m p o r t a n t p h y s i c a l l e n g t h scale of t h e p r o b l e m . T h e t i m e step At t h e n f o l l o w s b y considering a Courant-Friedrichs-Lewy ( C F L ) c o n d i t i o n : cAt < Aa;, or, f o r r e s t r i c t e d d e p t h ,

C F L = < 1 . (3.19)

T h e value of At depends o n t h e i n t e g r a t i o n m e t h o d . W e n o w take

E u l e r C F L < 0.5

A B M C F L < 0.8 .

W e use A B M i n t e g r a t i o n m o r e or less standardly.

For problems o f s o l i t a r y waves, t h e g r i d size m a y be determined o n t h e basis of t h e water d e p t h h. A n i n d i c a t i v e value of Ax/h f o r m o d e r a t e l y h i g h solitary waves is 1/5.

Hamiitonian formuiation H782 May 1994

T h e step-size m a y also be based o n the width of t h e s o l i t a r y wave. W i t h

'x - cV

C{x,t) = Hsech'

A (3.20) t h e w i d t h A can be e s t i m a t e d f r o m ^ . A (3.21)

As an example, f o r t h e case of a s o l i t a r y wave w i t h h e i g h t H = A m progressing o n w a t e r of d e p t h / i = 10 m , we have A = 18.3 m . Choosing Aa; = 2 m , we t h e n have 9 mesh p o i n t s over the w i d t h of t h e s o l i t a r y wave. Correspondingly, A t < 0.16 s.

For p e r i o d i c waves t h e smallest wave l e n g t h of interest should be considered. Consider, f o r example, a submerged bar w i t h the water d e p t h changing f r o m 0.40 m at the toe t o 0.10 m at the t o p and an i n c i d e n t t r a i n o f p e r i o d 2.02 s. U s u a l l y t h e t h i r d -h a r m o n i c component, w -h i c -h is generated t -h r o u g -h n o n l i n e a r i t y , is of i m p o r t a n c e f o r such problems. W e have t h e f o l l o w i n g wave lengths f o r t h e f o u r harmonics w i t h p e r i o d Ti = T / i at 0.10 m a n d 0.40 m d e p t h .

Ti

2.02 1.01 0.67 0.505

h

0.40 3.74 1.49 0.70 0.39 0.10 1.97 0.93 0.56 0.37

Table 3.1: Wave lengths Aj for the harmonics Ti = T/i, 1 , . . . , 4 .

W i t h Aa; = 0.10 m we have 5 p o i n t s per wave l e n g t h f o r t h e t h i r d h a r m o n i c at t h e shallower p a r t . A s e n s i t i v i t y test shows the difference between choices o f A a; = 0.1 m a n d Aa; = 0.05 m . For t h e smaller Aa; the representation of the smaller features seems t o be m o r e accurate t h a n f o r the larger Aa;. T h e general p i c t u r e of the t w o wave profiles is t h e same i n b o t h cases.

•^Foi solitary waves a.s resulting from different one-way equations, tlie precise form of tlie widtli A and the celerity c may differ, but these differences are not large (Dingemans, 1994b, section 6.3)

4 Spectral approach

I n t h e present chapter, we discuss spectral or pseudo-spectral s o l u t i o n o f t h e nonlinear water wave p r o p a g a t i o n based o n the explicit H a m i l t o n i a n m o d e l . One of the basic ideas b e h i n d a s p e c t r a l s o l u t i o n is t o i m p r o v e the n u m e r i c a l efiiciency of the s o l u t i o n b y achieving higher spatial r e s o l u t i o n w i t h a f e w spectral parameters. We h a d seen t h a t using a Fourier series over a f i n i t e span was n o t feasible i n connection w i t h the present H a m i l t o n i a n f o r m u l a t i o n s . O n t h e other h a n d , the constraint o f f i n i t e energy and the r e s t r i c t i o n of t h e i n i t i a l l y specified disturbance t o f i n i t e span p o i n t s t o t h e sinc-series as a n a t u r a l choice. W e discuss several possibilities of the n u m e r i c a l m e t h o d using the sine a p p r o x i m a t i o n .

4.1 Method A: a;-space

We recall f r o m section 3.1 t h a t the (regularised) e v o l u t i o n equations i n c l u d i n g the short-wave n o n l i n e a r i t y s are: 1 d \ Ci = - - T^[v(l + e)v]ix)+ / (ia;'A(a;',a;)logtanh|/(a;,3;') TT ax I J-oo (4.1a) - B l dx _J—oo r dx'[v{l^-e)nCt]^Bi[v{l-^e)7jCt]. (4.1b)

I n t h e simplest a d o p t i o n of sinc-series, C,{x) and v{x) i n (4.1a) a n d (4.1b) can be a p p r o x i m a t e d t h r o u g h the c a r d i n a l series:

i = o

CM = E C i W « i i ^ c ( ^ ) , (4.2a)

vix,t) = E ^ i W ^ i ^ M ^ l • (4-2b) i = o

F u r t h e r s i m p l i f i c a t i o n is achieved b y recognising t h a t the f u n c t i o n

A(a;', x ) l o g t a n h \f(x,x')\ belongs t o the PaleyWiener class of B(Ax) ( L u n d a n d B o w -ers, 1992; chapter 2) and hence can be expressed b y the sinc-series, i.e.

A(a;', a;) l o g t a n h |/(a;,a;')| = E A | ( a ; ) s i n c f ^ ^ ^ ^ ) (4.3) i = o V Aa; y

Hamiitonian formulation H782 May 1994

\*Ax)

=

0, [x = jAx]

( i ) - % X ^ * ^ K ^ ) V o g t a i i l i | / ( a ; , a ; , - ) | [ x ^ j A x ] . ( 4 . 4 )

[1 + S j j t j j I

I n t i o d u c i n g t h e expression (4.3) i n t o (4.1a) a n d using t h e i n t e g r a l r e l a t i o n {e.g., L u n d a n d Bowers, 1992; p 2 5 )

/•°° . f x ' - X j \ , ,

/ sine •—;—- dx = Ax

J-oo \ Ax J

yields the f o l l o w i n g s i m p l i f i e d e v o l u t i o n e q u a t i o n f o r surface elevation:

(4.5)

TT dx -7r[v{l+e)7]]{x) + AxJ2^j 3=0

(4.6)

T h e replacement of (4.1a) t h r o u g h (4.6) is i n f a c t t h e same as t h e t r a p e z o i d a l i n t e g r a -t i o n . However, -t h e f o r m a l i s m o f -t h e sine series under -t h e a s s u m p -t i o n -t h a -t -t h e f u n c -t i o n X{x', a ; ) l o g t a n h \ f { x , x')\ belongs t o t h e P a l e y - W i e n e r class o f B{Ax) gives j u s t i f i c a t i o n t o t h e exactness o f (4.6).

I n t h e i n t e g r a l over (—oo, x] i n (4.1b) we n o t e t h a t the c o m b i n e d f u n c t i o n v{l + e)r]Ci can b e a p p r o x i m a t e d t h r o u g h t h e o r t h o g o n a l sine series since v,(i E B{Ax) w h i l e ( 1 + e) a n d 77 are b o u n d e d ( 0 ( 1 ) ) . I n t r o d u c i n g t h e series

v{l + s)r)Ct{x') = Y^[v{l + s)vCi]j sine

i = o Aa; (4.7) i n t o ( 4 . 1 b ) results i n t h e f o l l o w i n g s i m p l i f i e d e v o l u t i o n e q u a t i o n f o r v: -gC. + Bl [v{l + e)r)Ct]

+

- B l dx

jy{l

+ eU,]i

r

dx' j=o smc X — x^ Ax (4.8) N o t e t h a t f o r x = A;Aa; a n d Xj = j ' A a ; t h e i n t e g r a l / dx smc — - — - = dx smc — — - -j- dx smc J-00 \ Ax J J-00 \ Ax J Jxi Ax Ax Ax r^^-j) sin 6» — + ^ r - ' ^ '-^dO ( 4 . 9 ) TT Jo

o

is t i m e i n v a r i a n t a n d can be easily c o m p u t e d . delft hydraulics 2 1

4.2 Method A: x - s p a c e

One o f t h e disadvantages of the m e t h o d using (4.6) is associated w i t h t h e e v a l u a t i o n of t h e i n t e g r a l f u n c t i o n f { x , X j ) appearing i n X j . T h i s can be avoided b y using t h e e v o l u t i o n equations i n x-space (Radder, 1992) given b y :

1 f°° TT

m = - dx'v(x')logtanli-\x-x'\ (4.10a)

TT J - o o 4

U n d e r the assumption t h a t t h e short-wave n o n l i n e a r i t y £ = 0, we have

1 r/r

Bl = ^ , .nd - = , . (4.11)

Expressing v[x) i i i sinc-series i n t h e u n i f o r m x - g r i d , the e v o l u t i o n equations become

m = E ^ i ^ i ( x ) (4-12a) i=o ^ = - ö C - ^ E E ^ > f c ^ i , A - ( x ) (4.12b) j=ok=o where Vj is t h e value of v at Xj = j^X a n d ^ i ( x ) = ^^X'sine l o g t a n h I | x - x ' | , (4.13a) A I smc

A x y V A x

(4.13b)

T h e integrals I j { x ) a n d I j j . { x ) are time-independent and can be accurately evaluated. N o t e , however, t h a t t h e e v o l u t i o n rates ih and cp are the E u l e r i a n e v o l u t i o n rates a n d do n o t therefore correspond t o the e v o l u t i o n rates at a fixed x - g r i d p o i n t . T h i s means t h a t a correspondence between the u n i f o r m x - g r i d and the s - g r i d has t o be made at each t i m e step unless the e v o l u t i o n rates of m and ip are derived at fixed x - g r i d p o i n t s .

4 . 2 . 1 E v o l u t i o n r a t e a t f i x e d x-gt'ids

Hamiltonian formulation H732 May 1994

dt ^ dt ^ dx dt

(4.14)

w h e r e dip/dt is t h e E u l e r i a n e v o l u t i o n r a t e o f ip (mass m or p o t e n t i a l (p). T h e convective t e r m can be expressed e n t i r e l y i n t h e x-space, i.e.;

a n d

dip _ dip dx _ dip dx dx dx Tj dx

dx] d dx , d r^i

dt dt J-oo dx

m dx, d [XJ , f^i . ,

/ -i-dx = T I / Vdx= Cidx •

J-oo ClX dt J-oo J-oo

(4.15a)

(4.15b)

4.3 Method B

I n t h i s p r o c e d u r e , one begins w i t h t h e expression f o r t h e H a m i l t o n i a n

/

oo 2^ ]^ roo roo rx

^gCdx - — dx J^^ dx'v(x)v(x')logt<inh- ^ dr

(4.16)

U s i n g the same sine series (4.2a)-(4.2b), t h e H a m i l t o n i a n H can be expressed as

^ = E e l A x - ; ^ E E ^ i ^ ' ^ ^ O ' ' ^ ' ^ ) (4-17) where

/

CO poo dx / dx' -oo J—oo smc Ax smc X - Xk l o g t a n h •

L

Ax ^' dr

T h e basic e v o l u t i o n equations are t h e n o b t a i n e d f r o m

1

dn

1

dn

m„ =

-Ax dvp' Ax dCp'

Hence, f r o m (4.17) and (4.19) one has

(4.18)

(4.19)

(4.20a)

Basic t o t l i e efficiency of t h i s approach is a s m a r t e v a l u a t i o n o f t h e i n t e g r a l r e l a t e d t o t h e k i n e t i c energy u s i n g t h e p r o p e r t y of t h e sine a n d l o g t a n h f u n c t i o n s .

4.4 Method C

Radder (1993a, a p p e n d i x A ) has s h o w n t h a t T , t h e k i n e t i c energy p a r t o f the t o t a l H a m i l t o n i a n , can be expressed i n t h e f o r m :

r =

Y.Y.^>mj-k\-Ax)

(4.21) j=0k=0

where

a n d Vj denotes v at t h e j-th p o i n t along t h e equidistant x g r i d or i n other words,

<x) = T . ^ ] s m c ( ^ ^ y (4.22b)

N o t e t h a t — k\; Ax) is exact (evaluated e x a c t l y on the x-space w i t h o u t using any a p p r o x i m a t e t r a n s f o r m a t i o n t o a;-space) a n d needs t o be evaluated o n t h e x - g r i d once at the s t a r t o f t h e c o m p u t a t i o n . T h e disadvantages of u s i n g (4.22a), however, become apparent i f t h e e v o l u t i o n equations are l o o k e d at.

4 . 4 . 1 E v o l u t i o n e q u a t i o n : 7]

Since t h e pseudo-spectral variables v^'s are n o t canonical, t h e e v o l u t i o n rates cannot be expressed as s i m p l e f u n c t i o n a l derivatives of the H a m i l t o n i a n i n t e r m s o f Vj's; i.e., t h e e q u a t i o n

^* = - ^ . j i = - 2 E ^ K I i - ^ ' l ; A x ) (4.23) _V

Ax Sv;

i s n o t v a l i d . O n t h e other h a n d , t h e e v o l u t i o n e q u a t i o n f o r mass (consequently, surface elevation) s h o u l d be d e t e r m i n e d f r o m

. L \ - ^ Z 1 M (A9A^

AxSvp- Ax^Sv^Svj,- ^

Hamiltonian formulation H782 May 1994

k \ J

wliere rcj- denotes t i e ^-coordinate of the j - t h g r i d p o i n t o n the x-space.

4 . 4 . 2 E v o l u t i o n e q u a t i o n : cp

T h e expression

is exact f o r the p o t e n t i a l energy V where ( j denotes the surface elevation f r o m the s t i l l water level at t h e j-th p o i n t along the u n i f o r m a ; - g r i d . Using (4.26) and (4.21) the e v o l u t i o n e q u a t i o n f o r ipp becomes

f p = -sCv + ^ E E ^ d i - k\;Ax)^^v*vl . (4.27)

We have t r e a t e d the i n t e g r a l I{\j - k\; Ax) i n (4.27) i n v a r i a n t (on u n i f o r m

x-gi'id)

leaving the discrete t»|'s t o respond t o increments i n surface elevation ( j ( o n u n i f o r m x - g r i d ) . T h e derivative d{vjVl)/dCp can be c o m p u t e d i n the f o l l o w i n g way.

L e t xtj denote the j—th p o i n t on a n o n n n i f o r m l y d i s t r i b u t e d x - g r i d such t h a t a;*- = x ( x j ) a n d Aa;*- represent the s h i f t i n a;*- due t o A ^ p . Avj(Cp), the i n c r e m e n t i n v*j due t o ACp, is t h e n given b y

Av^iCp) = v(x) + Ax}) - vix}) = ^ x } ^ i x } ) +

T h u s ,

TrrV-iVi, = l i m dCp ^ ^ ACp^o

T h e order t o w h i c h d e p t h v a r i a t i o n and surface n o n l i n e a r i t y are e x p l i c i t l y accounted f o r i n the m o d e l depends on the f o r m of the a p p r o x i m a t i o n of the a c t u a l Woods t r a n s f o r m a t i o n . I n f i r s t instance, we shall use the r e l a t i o n = t] w h i c h gives

dx x{x) — x{xq) + / rjdx- A n i n t e g r a l expression f o r dx^/d(p t h e n results: Jx d x \ . . dri , r ) d r ] l , P i 1 . f x - X p \ ^ ,^ 1 ^ ^ ) = / 1^'^X= / •7~-dx= / - s m c — - ^ ] d x . (4.28) dCp Jx dCp Jxl dCpV 7? V Aa; y delft hydraulics ₂₅

4.5 Summary

T l i e procedure o u t l i n e d i n M e t h o d A eliminates a significant a m o u n t of n u m e r i c a l c o m p u t a t i o n b y u t i l i s i n g the properties of sine f u n c t i o n s . Besides, one expects t h a t smaller number of g r i d p o i n t s are r e q u i r e d t h a n t h a t i n the previous t i m e - d o m a i n m o d e l t o achieve the same degree of accuracy. T h e n u m e r i c a l efficiency o f M e t h o d B depends largely on how t h e i n t e g r a l E{j,k,'q) can be f u r t h e r s i m p l i f i e d .

T h e pxocedure used i n M e t h o d C has one d i s t i n c t advantage f r o m the physical p o i n t of view: the e v o l u t i o n equations are kept exact, the a p p r o x i m a t i o n being i n t r o d u c e d o n l y t h r o u g h the t r a n s f o r m a t i o n f r o m a; t o x a n d vice versa. T h i s allows easier m o d i f i c a t i o n of t h e n u m e r i c a l code t o i n c o r p o r a t e higher-order f o r m u l a t i o n s . T h e disadvantage is r e l a t e d t o t h e f a c t t h a t the n u m e r i c a l operations are d i s t r i b u t e d over u n i f o r m x - g r i d and x-grid i n v o l v i n g continuous i n t e r p o l a t i o n f r o m one t o another.

Hamiitonian formuiation H782 IVIay 1994

5 Numerical method using sinc-series

I n t h e previons chapter we discussed possible f o r m u l a t i o n s based on sine m e t h o d of w h i c h we set out the reasons f o r the choice of the e v o l u t i o n equations i n the x-space of M e t h o d A (repeated here f o r convenience):

Ct = 1 d TT dx -ir[v{l+e)r]]{x) + AxY^X] 3=0 Vt = -gCx + Bi[v{l + e)r]Ct] ^ " ^ E K l + e)vCt]j r dx' • n J — OO

+

_dx 3=0 smc , oo V A x (5.1) (5.2)

w i t h A* defined i n (4.4). N o t e t h a t t h e choice of the sinc-series f o r t h e g l o b a l ap-p r o x i m a t i o n requires t h e g r i d ap-points t o be u n i f o r m l y sap-paced a n d f o r x^ = kAx and Xj = jAx, one has

, , . f x ' - x , - \ A x A x r / dx smc — - — - = —— H / y_oo

V A x y 2 TT Jo

Ax . A x

ri''-^)

s i n ö

1"

de.

(5.3)

5.1 Numerical procedure

T h e n u m e r i c a l procedure proceeds i n t w o steps: c o m p u t a t i o n of the e v o l u t i o n rates a n d t h e t i m e i n t e g r a t i o n of t h e e v o l u t i o n . T h e r e are t w o available options i n t h e code f o r c a r r y i n g out t h e t i m e i n t e g r a t i o n : a f i r s t order e x p l i c i t m e t h o d ( E u l e r ) a n d a f o u r t h order predictivecorrective m e t h o d ( A B M ) . I t has been seen f r o m previous c o m p u t a -tions using t h e t i m e - d o m a i n m o d e l t h a t m u c h larger t i m e steps can be successfully used w h i l e using the A B M m e t h o d f o r t h e t i m e i n t e g r a t i o n of t h e e v o l u t i o n rates. T h e t i m e - i n t e g r a t i o n procedure is the same i n b o t h t h e models.

Before describing some specific aspects of t h e c o m p u t a t i o n of the e v o l u t i o n r a t e of sur-face elevation, we l o o k at t w o m a i n n u m e r i c a l operations, i n v o l v e d i n t h e c o m p u t a t i o n of t h e e v o l u t i o n equations. These are i n t e g r a t i o n of a f u n c t i o n over an i n t e r v a l a n d t h e e s t i m a t i o n o f t h e derivative at a n o d a l p o i n t .

5 . 1 . 1 I n t e g r a t i o n p r o c e d u r e

Consider the evaluation of

Ig = g{x)dx (5.4)

Jxi

where Xi is t h e x - c o o r d i n a t e o f t h e g r i d p o i n t i and g is an a r b i t r a r y regular f u n c t i o n . A Gauss-quadrature procedure is i m p l e m e n t e d i n order t h a t t h i s i n t e g r a l can be evaluated accurately f o r r e l a t i v e l y large êx. T h i s means t h a t t h e f u n c t i o n value needs t o be k n o w n at the i n t e g r a t i o n p o i n t s i n t h e i n t e r v a l [xi, S j + i ] . E s t i m a t i o n of these values is done t h r o u g h a higher-order a p p r o x i m a t i o n based on shape-function m e t h o d (see Zienkiewicz and M o r g a n , 1983 f o r the d e f i n i t i o n of shape f u n c t i o n s ) . I n t e r m s of the shape f u n c t i o n s A^;(/i)'s defined i n parameter /x, an a p p r o x i m a t i o n t o the f u n c t i o n g{x) is expressed i n t h e f o r m l=2m ^ ( « ( M ) ) = E ff(«(i-m+o)^K/^) (5.5) One can n o w r e w r i t e t h e i n t e g r a l Ig as

/=2m fl-tixi+i) dx

h=2Z9{i-m+l)l^ Ni{p)-r-dii. (5.6)

N o t e t h a t the order of t h e a p p r o x i m a t i o n is ( 2 m — 1) a n d t h e a p p r o x i m a t i o n of t h e f u n c t i o n over Xi < x < x^^i is based o n the discrete values of the f u n c t i o n at nodes s u r r o u n d i n g the i n t e r v a l ( i . e , , m number of nodes on either side of t h e i n t e r v a l are used f o r the a p p r o x i m a t i o n ) . I n a general way f o r b o t h u n i f o r m and n o n u n i f o r m grids m a p p i n g of a; t o yti can t a k e place i n a similar way ( k n o w n as 'isoparametric' approach) i n t h e f o r m

X{fi) = E ^{i-m+l)Nl{l^) • (5.7)

1=1

A simpler t r a n s f o r m a t i o n v a l i d o n l y f o r u n i f o r m l y spaced grids is

a;(/«) = + ^ ^ {fx + 1) . (5.8)

I t is convenient t o express t h e i n t e g r a l (5.6) over the i n t e r v a l [-1,1] f o r t h e Gauss-qua-d r a t u r e integrations. T h i s is Gauss-qua-done b y another t r a n s f o r m a t i o n , Gauss-qua-defineGauss-qua-d b y

^ = ^ ( , , ) + / < £ i ± l t i M ( f + i ) , ( 5 . 9 ) where M(xO = - ^ , M ( ^ ( W = ^ - (5-10) F i n a l l y , one has '=2™

nl dx da

= E Hi-m+l) ]_^^^-d{xW^^ Si-m+lll (5.11)

where

Hamiltonian formulation H782 May 1994

(5.12)

w h i c h is c o m p u t e d t h r o u g h Gauss-quadrature, is c o m p u t e d o n l y once at the s t a r t of t h e c o m p u t a t i o n .

5.1.2 C o m p u t a t i o n o f t h e d e r i v a t i v e

T h e derivative o f a f u n c t i o n at a g r i d p o i n t i is c o m p u t e d b y i n t e r p o l a t i o n f r o m i t s values at the adjacent g r i d p o i n t s f r o m the r e l a t i o n :

J2 wjf{i + l) (5.13)

where the coefficients w;'s are independent of i ( t h o u g h dependent o n Aa;) f o r u n i f o r m grids. T h e order of i n t e r p o l a t i o n is 2m (based on 2m + 1 nodes).

5 . 1 . 3 E v o l u t i o n rates

C o m p u t a t i o n of t h e e v o l u t i o n r a t e of v f r o m (5.2) is r a t h e r s t r a i g h t - f o r w a r d using the i n t e g r a t i o n and derivative procedures j u s t described. C o m p u t a t i o n o f the e v o l u t i o n r a t e o f surface elevation needs some specific considerations.

C o m p u t a t i o n of t h e t e r m A | (see E q . (4.4)) involves the expression

l o g t a n h | / ( a ; f c , a ; j ) | = l o g t a n h

L

^* dx ( l + £)7? • T h e i n t e g r a l argument of l o g t a n h is evaluated t h r o u g h where / = 1,2,. (5.14) (5.15) _ dx , n (5.16)

I t is f u r t h e r t r u e t h a t the value of decreases r a p i d l y away f r o m t h e p o i n t i due t o t h e l o g t a n h operator. I t is possible therefore t o l i m i t t h e s u m t o smaller n u m b e r of nodes, e.g.,

3=n

j=0 3=i—p

(5.17)

0.

-2 -1. 0. 1. 2. Figure 5.1; The function logtanh(a;)

where the m i n i m u m value of (i - p) a n d the m a x i m u m value of (i + p) are l i m i t e d b y the c o m p u t a t i o n a l g r i d p o i n t s used. T h e i n d e x ^) at a g r i d - p o i n t i is calculated b y using a simple c r i t e r i o n based o n the l o c a l water d e p t h

pAx

> Rl (5.18)

where Ri is an user specified l i m i t . T h e value of Ri should never be lower t h a n 3. I t i s , however, safer t o use a somewhat larger value (say a r o u n d 10). T h e behaviour of the f u n c t i o n logtanh(a;) is shown i n F i g . 5 . 1 , where i t has t o be u n d e r s t o o d t h a t at a; = 0 an a s y m p t o t e is present. Here t h e g r a p h is such t h a t we stay clear f r o m a; = 0 b y 0.005 u n i t s .

5 . 1 . 4 E n d t r e a t m e n t

T h e basic f o r m u l a t i o n s w h i c h c o n s t i t u t e t h e n u m e r i c a l m o d e l here are based o n an i n f i n i t e span. T h e v a l i d i t y of a finite c o m p u t a t i o n a l l e n g t h rests o n the f a c t t h a t t h e entire wave a c t i o n is contained w e l l w i t h i n the finite span. T h i s is difRcult t o meet i n practice i f the c o m p u t a t i o n s are t o be r u n f o r a l o n g t i m e of wave p r o p a g a t i o n . A proper p h y s i c a l t r e a t m e n t o f t h i s p r o b l e m , t h o u g h r a t h e r i m p o r t a n t , is n o t y e t u n d e r t a k e n . I n t h e f o l l o w i n g , n u m e r i c a l aspects of t h e 'end t r e a t m e n t ' are b r i e f l y described.

T h e procedures used f o r c o m p u t i n g the integrals and derivatives as described earlier assume t h a t t h e required number of g r i d p o i n t s are available o n either side ( o f an i n t e r v a l or a g r i d - p o i n t ) and are therefore n o t s t r i c t l y v a l i d at t h e ends. T h i s means t h a t the procedures need t o be m o d i f i e d at the ends. A n a t t r a c t i v e a l t e r n a t i v e is t o define the necessary n u m b e r of a r t i f i c i a l g r i d - p o i n t s t o be able t o use the same procedures f o r a l l t h e r e a l nodes. Presently, the f u n c t i o n value ƒ at an a r t i f i c i a l g r i d -p o i n t x"- is defined t h r o u g h the sim-ple r e l a t i o n :

/(^") = f i x ' ) (5.19)

where x'^ is t h e closest real g r i d - p o i n t . T h e a r t i f i c i a l g r i d - p o i n t s are also assumed t o be equidistant i n t h e same manner as the r e a l g r i d p o i n t s . T h o u g h the simple f o r m (5.19) can be m o d i f i e d , i t is n o t considereded r e a l l y i m p o r t a n t w i t h t h e a s s u m p t i o n t h a t the wave a c t i o n should be non-significant near the ends.