Numerical modeling of Boussinesq equations by finite element method

Coastal Engineering 37 Ž1999. 97–122 www.elsevier.comrlocatercoastaleng

Numerical modeling of Boussinesq equations by finite element method Y.S. Li a

a,)

, S.-X. Liu a , Y.-X. Yu b, G.-Z. Lai

b

Department of CiÕil and Structural Engineering, The Hong Kong Polytechnic UniÕersity, Hong Kong, China b State Key Laboratory of Coastal and Offshore Engineering, Dalian UniÕersity of Technology, Dalian 116024, China Received 15 June 1998; received in revised form 27 January 1999; accepted 4 February 1999

Abstract In this paper, a numerical model based on the improved Boussinesq equations derived by Beji and Nadaoka wBeji, S., Nadaoka, K., 1996. A formal derivation and numerical modeling of the improved Boussinesq equations for varying depth. Ocean Eng. 23 Ž8., 691–704x is presented. The finite element method is used to discretize the spatial derivatives. Quadrilateral elements with linear interpolating functions are employed for the two horizontal velocity components and the water surface elevation. The time integration is performed using the Adams–Bashforth–Moulton predictor–corrector method. Five test cases for which either theoretical solutions or laboratory results are available are employed to test the proposed scheme. The model is capable of giving satisfactory predictions in all cases. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Boussinesq equations; Finite element; Absorbing boundaries

1. Introduction The ambient wave conditions must be determined before the design of coastal structures can proceed. An accurate estimation is essential for human safety and cost-effective design. As offshore waves propagate inshore into shallow water, they are transformed as a result of the variation in the bottom topography and the presence of barriers such as islands and breakwaters. Because of non-linear effects, there are energy transfers among the different components of shoaling waves. Realistic modeling of the transformation process should include )

Corresponding author. Tel.: q852-2766-6069; Fax: q852-2334-6389

0378-3839r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 3 9 Ž 9 9 . 0 0 0 1 4 - 9

98

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

these non-linear interactions. A review of the wave transformation in the nearshore zone was given by Hamm et al. Ž1993.. At present, the Boussinesq-type equations, which are weakly non-linear and dispersive, have been found to give a relatively accurate description of wave transformation in coastal regions, if the underlying assumptions in their derivation are not violated. These depth-integrated equations result from the simplification of a three-dimensional problem to an equivalent two-dimensional one, making their solution possible because of significant savings in computation time. The first such set of equations for variable water depth was derived by Peregrine Ž1967., which are referred to as the standard Boussinesq equations and became a starting point for many studies of the same kind. Unfortunately, errors in the prediction of waves occur when the standard Boussinesq equations are applied in intermediate or deep water, primarily due to the inaccuracy of the linear dispersion relation with increasing water depth. To extend their range of applicability, efforts were made by a number of investigators to improve the standard equations or derive alternative forms. Based on the Pade´ approximant of the exact linear dispersion relation, Witting Ž1984. introduced a formulation with improved dispersion characteristics and gave a new set of Boussinesqtype equations for a single horizontal dimension. Madsen et al. Ž1991. presented a new form of equations in terms of the depth-integrated velocities for two horizontal dimensions with improved linear shoaling and dispersive characteristics and Madsen and Sørensen Ž1992. rederived the new Boussinesq equations for slowing varying bathymetry. Beji and Nadaoka Ž1996. gave a similar set of Boussinesq equations with exact energy conservation by adding and subtracting a dispersive term in the momentum equations. Nwogu Ž1993. derived alternative sets of Boussinesq equations in terms of the velocities at an arbitrary level. All these equations have similar dispersive properties with their corresponding selected values of the parameters, but different shoaling characteristics. In particular, the shoaling characteristics of the equations given by Madsen and Sørensen Ž1992. is acceptable up to kh s 3 ŽSchaffer and Madsen, 1998.. A detailed review of the ¨ wave transformation formulations over uneven bottoms was recently given by Dingemans Ž1997.. Parallel to the formulations in the time domain described above, frequency domain formulations have been developed also. Examples are from Kirby Ž1990. for the frequency domain version of the conventional equations, Mase and Kirby Ž1992. and Madsen and Sørensen Ž1993. for extended equations with improved linear shoaling and dispersive characteristics. Up to now, most of the numerical models based on Boussinesq-type equations in the time domain employ the finite difference method ŽAbbott et al., 1984; Wei and Kirby, 1995; Beji and Nadaoka, 1996.. The finite difference method is easy to use, but is not versatile enough to deal with irregular boundaries which occur in many coastal engineering problems. In contrast to the finite difference method, the finite element method, though it requires more programming efforts, can handle geometrically complex domains because it is based on unstructured grids and has been applied broadly in coastal engineering. Another advantage of the finite element method is the possibility to minimize the number of grid points by using only a fine resolution in the areas of interest. For example, the finite element method has been used to solve the non-linear shallow water equations ŽKashiyama and Ito, 1995; Li and Zhan, 1998; etc... However,

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

99

the use of the finite element method to solve Boussinesq-type equations is rare. Antunes do Carmo et al. Ž1993. and Antunes do Carmo and Seabra Santos Ž1996. solved the Boussinesq equations using the finite element method, excluding the dispersion modification terms. In this paper, a numerical model using finite element method, based on the Boussinesq-type equations given by Beji and Nadaoka Ž1996., is presented. The third-order spatial derivatives of the water surface elevation that appear in the dispersion modification terms require special treatment. This will be described in detail in Section 3. Solitary wave propagation for a long distance is simulated to test the stability and conservative properties of the model. Results from the computation of the wave focusing behind a topographical lens are compared with the laboratory data given by Whalin Ž1971. to investigate the generation of higher harmonics by varying water depth. Finally, the wave diffraction through a breakwater gap is modeled and the results are compared with the experimental results of Yu et al. Ž1998.. In Section 2, the governing equations are briefly described. The finite element scheme is presented in Section 3 and five test cases are reported in Section 4. Lastly, the conclusion is given in Section 5.

2. Governing equations The Boussinesq equations in terms of the depth-averaged velocity with improved dispersive characteristics derived by Beji and Nadaoka Ž1996. are as follows:

ht q = P Ž h q h . u s 0

Ž 1.

u t q Ž u P = . u q g=h s Ž 1 q b .

h 2

= = P Ž hu t . q b

yŽ 1 q b .

h2 6

gh 2

= = P Ž ut . y b

= = P Ž h=h . gh2 6

= = 2h

Ž 2.

where u s Ž u,Õ . is the two-dimensional depth-averaged velocity vector, h is the water surface elevation, h s hŽ x, y . is the water depth as measured from the still water level, g is the gravitational acceleration and b is a constant to be given below. The subscript t denotes partial differentiation with respect to time. The dispersion relation of Eqs. Ž1. and Ž2. after linearization for uniform water depth is: c2 s

v2 k2

s

gh Ž 1 q b k 2 h2r3 . 1 q Ž 1 q b . k 2 h2r3

Ž 3.

where c is the wave celerity; v is the wave circular frequency and k is the wave number. The constant b can be determined by matching Eq. Ž3. with the second-order Pade´ approximant of the full linear dispersion relation v 2rgk s tanh kh given below.

v2 s gk

kh Ž 1 q k 2 h 2r15 . 1 q 2 k 2 h 2r5

.

Ž 4.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

100

Fig. 1. Comparison of dispersion relation for different values of b with linear theory.

Comparing Eq. Ž3. with Eq. Ž4., it is seen that b s 1r5 would make the two equations identical. In fact, b can take other values in different asymptotic expansions of the linear dispersion relation. Fig. 1 is the comparison of the dispersion relation, Eq. Ž3., for different values of b with the linear theory. It can be seen that the one corresponding to the Pade´ approximant is the best alternative. So, b s 1r5 is used in this study. This choice corresponds to the improved Boussinesq equations of Madsen and Sørensen Ž1992. with b s 1r15 and those of Nwogu Ž1993. with a s y2r5 and extends the applicability of the Boussinesq equations to relatively deeper water. Eqs. Ž1. and Ž2. are used in this paper because of exact energy conservation. However, it should be noted that the shoaling of the equations is only accurate up to kh s 1.7.

3. Numerical model 3.1. Finite element formulation Grouping the time derivatives together, Eqs. Ž1. and Ž2. can be re-written in scalar form as follows: E

ht q

Ex

E

Ž Hu . q

Eu rt q u

Ex

6

Ž HÕ . s 0 Eh

Eu qÕ

b gh2 q

Ey

ž

Ey

qg

E2 c Ex2

Ž 5.

b gh E y

Ex E2 d

q E xE y

2

/

s0

Ex

E Ž hc . Ex

q

E Ž hd . Ey

Ž 6.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

EÕ st q u

Ex

b gh2 6

Õ y b1

Ey

qg

Ey

q E Ž hu .

2 Ex

Ex

h E

E Ž hu .

2 Ey

E Ž hc .

E y2

/

Ey E Ž hÕ .

q

Ex

q

E Ž hd .

Ex

Ey

s0

E Ž hÕ .

q

Ey

2

E2 d

E xE y

h E

b gh E y

E2 c

ž

q

u y b1

Eh

EÕ qÕ

Ey

101

Ž 7. q b1

q b1

h2 6 h2 6

ž ž

E2 u

E2 Õ q

Ex2

E xE y

E2 u

E2 Õ q

E xE y

E y2

/ /

sr

Ž 8.

ss

Ž 9.

in which b 1 s 1 q b , H s h q h , c s EhrE x and d s EhrE y. Eqs. Ž6. – Ž9. can be simplified by neglecting the second-order spatial derivatives of h which are small for slowly varying water depth. Eu rt q u

Ex

Eh

Eu qÕ

Ey

qg

Eh Ec y b gh

y Ex Ex

EÕ st q u

Ex

Ey

b1 3

h2

Eh

y Ey Ey E2 u Ex

2

3

y

E2 c Ex2

bg y 3

h

2

E2 d E xE y

bg y

Ey

3

h2

E2 c

bg y

E xE y

3

h2

h2

E2 Õ E xE y

Eh Eu y b1 h

y Ex Ex

Ž 10 .

E2 d E y2

b g Eh Ec b g Eh Ec h y h s0 2 Ex Ey 2 Ey Ex

b1 3

h

2

b g Eh Ed b g Eh Ed h y h s0 2 Ey Ex 2 Ex Ey

qg

Eh Ed y b gh

uy

y Ex

EÕ qÕ

bg

Ž 11 .

b 1 Eh EÕ b 1 Eh EÕ h y h sr 2 Ey Ex 2 Ex Ey

Ž 12 . Õy

b1 3

h2

E2 u y E xE y

b1 3

h2

E2 Õ Ey

2

Eh EÕ y b1 h

y Ey Ey

b 1 Eh Eu b 1 Eh Eu h y h s s. 2 Ex Ey 2 Ey Ex

Ž 13 . Subdividing the domain of interest V into finite elements and denoting a generic element by V e , the dependent variables in Eqs. Ž5., Ž10. – Ž13. are approximated within the elements as follows: nd

f f Ý Ni f i

Ž 14 .

is1

where the f i are the values of any independent variable at the nodal points, nd is the number of nodes and Ni is the interpolation function.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

102

The Galerkin method is applied to Eqs. Ž5., Ž10. – Ž13.. As an example, the integral form of Eq. Ž5. can be written as:

HVHW

k

ž

E

ht q

E

Ž Hu . q

Ex

Ey

/

Ž HÕ . d xd y s 0

Ž 15 .

where Wk is the weighting function. Assuming Wk s Nk and substituting Eq. Ž14. for the variables into Eq. Ž15., the following matrix equation for an element V e can be obtained: Mk ih˙ i q Q k i j Hi u j q R k i j Hi Õj q Q k i j u i Hj q R k i j Õi Hj s 0

Ž 16 .

where the overdot denotes the partial differentiation with respect to time and: Mk i s

HV HN N d xd y ; k

e

Qk i j s

HV HN N k

e

Ž 17 .

i

ENj

i

Ex

d xd y ;

Rki j s

HV HN N k

e

ENj

i

Ey

d xd y.

Ž 18 .

Similar integral equations for Eqs. Ž10. – Ž13. can also be written. After some manipulations, including the use of integration by parts and Green’s theorem for the higher order derivatives, the following element matrix equations for an element V e can also be obtained: Mk i r˙i q Q k i j u i u j q R k i j Õi u j q Sk ihi q b g Ž A k i c i q Wk i d i . s 0

Ž 19 .

Mk i s˙i q Q k i j u i Õj q R k i j Õi Õj q Tk ihi q b g Ž X k i c i q Bk i d i . s 0

Ž 20 .

Ž Mk i q b 1 A k i . u i q b 1Wk i Õi s Mk i ri

Ž 21 .

Ž M k i q b 1 Bk i . Õ i q b 1 X k i u i s M k i s i

Ž 22 .

where: Sk i s g

ENi

HV HN

Ak i s Bk i s Wk i s

k

e

1

Ex

H Hh 3 V

2

H Hh 3 V

2

H Hh 3 V e

ENk ENi

2

ENk ENi

Ey Ex

H Hhh 2 V e

x

Nk

Ey

1

H Hhh 3 V

d xd y y

1

H Hhh 3 V e

ENk ENi

1 y

d xd y y

ENi

k

e

x

e

Ey Ey

e

1

HV HN

Ex Ex

e

1

d xd y ; Tk i s g

d xd y q

ENi Ey

1

d xd y y SWk i ;

Nk

y Nk

H Hhh 6 V e

d xd y ;

y Nk

ENi Ex ENi Ey ENi Ex

Ž 23 .

d xd y y SA k i ;

Ž 24 .

d xd y y SBk i ;

Ž 25 .

d xd y

Ž 26 .

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

Xk i s

1

H Hh 3 V

ENk ENi

2

Ex Ey

e

1 y

H Hhh 2 V e

y Nk

d xd y q

1

H Hhh 6 V e

ENi

x

Nk

ENi Ey

103

d xd y

d xd y y SX k i .

Ex

Ž 27 .

In the above equations, h x and h y are the derivatives of the water depth with respect to spatial coordinates x and y, respectively. SA k i , SBk i , SWk i and SX k i are boundary integrals given as follows: SA k i s SWk i s

1

ENi

2

Hh N 3 s

k

e

1

2

Hh N 3 s

k

e

Ex

1

n x d s;

SBk i s

n y d s;

SX k i s

ENi Ex

ENi

2

Hh N 3 s

k

e

1

2

Hh N 3 s e

k

Ey ENi Ey

n y d s;

Ž 28 .

n x d s.

Ž 29 .

where s e denotes the boundary of the element with unit normal n. These boundary integrals are needed only when natural boundary conditions are specified on part of the computational domain boundaries. In this model, linear quadrilateral elements are used. Since the first-order spatial derivatives are not continuous across adjacent elements, special care has to be taken. In this study, the nodal values of the first spatial derivatives of the surface elevation h at a particular node, c i and d i in Eqs. Ž19. and Ž20., are determined from the weighted average of the derivatives in the elements surrounding the node, i.e.: ci s

1

Nie

Ý

AeZ i js1

Aei j

Eh

ž / Ex

;

di s

ij

1

Nie

Ý Aei j

AeZ i js1

Eh

ž / Ey

Ž 30 . ij

where Nie is the number of elements surrounding the node i, Aei j the area of the jth e i element and AeZ i s Ý Njs1 Aei j the total area of the elements. ŽEhrE x . i j and ŽEhrE y . i j are the first-order spatial derivatives of the surface elevation with respect to x and y at node i calculated in the jth element, respectively. This corresponds to the central difference approximation of the first derivative of h if rectangular elements are used. 3.2. Predictor–corrector method The derivatives in time are discretized using the Adams–Bashforth–Moulton predictor–corrector scheme similar to that used by Wei and Kirby Ž1995.. To illustrate the method succinctly, the assembled global equations corresponding to Eqs. Ž16., Ž19. and Ž20. are written in matrix form.

w M x  f˙4 s  E 4

Ž 31 .

where f s h , r or s. M is the coefficient matrix and E is a known vector calculated from the known values of h , u, Õ, c and d. If the variable f is discretized at t s nD t, the predictor step is the third-order explicit Adams–Bashforth scheme.

w M x  f˜ 4

nq 1

n

s w M x f 4 q

Dt 12

n

23  E 4 y 16  E 4

ny1

q 5 E 4

ny2

Ž 32 .

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

104

where the values on the right hand side are known from calculations at previous time steps. The resulting set of simultaneous algebraic equations are denoted by: w M x  x4 s  F4 Ž 33 . where  x 4 is the unknown vector  f˜ 4 nq 1 , F is the constant vector corresponding to the right hand side of Eq. Ž32.. Eq. Ž33. are then solved using the following iteration procedure:

w M˜ x  x 4 lq1 s Ž w M˜ x y w M x .  x 4 l q  F 4 Ž 34 . in which l is the number of iteration and M˜ is a diagonal matrix. There are several ways of constructing diagonal mass matrices, known as lumped mass matrices ŽHughes, 1987.. Here, the concept of proportional lumping technique is used and the elements m ˜ ki of M˜ are determined as follows: m ˜ ki s0

Ž k/i. ;

Ž 35 .

m ˜ k k s a mk k where: ND ND

as

Ý

ND

Ý mk i

Ý mk k

ks1is1

Ž 36 .

ks1

in which ND is the total number of nodes. The iteration continues until the maximum difference in x i between consecutive iterations is less than 0.001. Usually, the number of iterations required is less than 10. The value of h˜ at time level Ž n q 1. can be obtained directly by solving Eq. Ž32.. However, the values of u˜ and Õ˜ have to be obtained by solving Eqs. Ž21. and Ž22. from the known values of r˜ and s˜ at time level Ž n q 1.. Because u˜ and Õ˜ are coupled, Eqs. Ž21. and Ž22. are solved by the following iteration procedure starting with the estimated values of Õ˜ at time level Ž n q 1. until convergence is achieved. Õ ™ Eq. Ž 21 . ™ u ™ Eq. Ž 22 . ™ Õ ™ Eq. Ž 21 . ™ u ™ . . . Ž 37 . The global equations in each iteration are solved using the bi-conjugate gradient ŽBICG. method ŽHoward et al., 1990.. The errors in u˜ and Õ˜ are defined as:

Ý Ž u˜ i y u˜Ui . D u˜ s

i

u 2i

ݘ i

2

Ý Ž Õ˜i y Õ˜iU . DÕ˜ s

i

ÝÕ˜i2

2

Ž 38 .

i

where u˜ i and Õ˜i are the current values and the U denotes the previous estimate. The iteration continues until the maximum value of D u˜ and DÕ˜ is less than 0.001. In general, only one iteration is needed. After the predicted values of h˜ , r, ˜ s,˜ u˜ and Õ˜ at time level Ž n q 1. are evaluated, the vector E˜ nq 1 can be calculated. Then the fourth-order Adams–Moulton scheme is used to calculate the final values of the variables. Dt nq 1 n n ny1 ny2 w M x  f 4 nq1 s w M x  f 4 q 9  E˜ 4 q 19  E 4 y 5  E 4 q  E4 24 Ž 39 .

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

105

The same procedure described in Eq. Ž37. is used to evaluate the final values of u and Õ at time level Ž n q 1.. In this model, all the algebraic equations are solved using iteration methods. To save computer storage which is a typical problem of the finite element method, only the non-zero elements of the two-dimensional matrices are stored. To ensure stability of the numerical scheme, the following criterion, which is the well-known CFL condition, must be satisfied: Dt cw F 1, Ž 40a. D lmin where c w is the wave velocity and D lmin is the minimum length of the sides of the finite elements. To ensure sufficient accuracy of the computed results, the following widely accepted criterion for wave simulation is applicable: L G 8 ; 10, Ž 40b. D lmax where L is the local wave length and D lmax is the maximum length of the sides of the finite elements. More than 10 grid points should be used per wave length for highly non-linear waves because of the generation of higher harmonics. The five test cases described in Section 4 verified the conditions Ž40a. and Ž40b.. In general, the computational cost of the finite element method is higher than that of the finite difference method, but this is compensated by its versatility in handling complex domains. For example, in the case of wave runup around a circular cylinder involving 4964 elements described in Section 4.5, the computer time required to advance the solution by one time step using a Pentium II Ž266 MHz. personal computer is 4 s. 3.3. Boundary conditions The conditions at the boundaries enclosing the computational domain must be specified to completely define the problem. Generally, there are four kinds of boundaries. 1. Incident wave boundaries which generate the waves of interest. 2. Fully reflective boundaries which reflect the arriving wave energy completely. 3. Absorbing wave boundaries or open boundaries which absorb all wave energy arriving from within the fluid domain. 4. Partially reflective boundaries which only reflect partly the arriving wave energy. The first three kinds of boundaries are considered below as they are required in Section 4 for the test cases. 3.3.1. Incident waÕe boundaries If the incident wave elevation h I is given, we can use linear theory to obtain the incident wave velocities as: v v u I s h I cos u Õ I s h I sin u Ž 41 . kh kh where u is the wave propagation direction relative to the x-axis.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

106

It should be noted that in the finite element scheme described above, the first-order derivatives of the surface elevation are calculated by Eq. Ž30.. If it is applied at the boundaries, only two surrounding elements are available for the calculation of the weighted average and this may lead to instabilities. In general, the gradient of the surface elevation should be specified at the incident wave boundary. Generally, the initial water surface elevation is set to zero for the whole domain. When the computation begins, the abrupt changes of the water surface at the incident boundary will cause some high-frequency oscillations of the wave surface elevation. To reduce these oscillations, the first two or three incident waves are divided by m defined by Eq. Ž44. below in which the parameter D d and d s now denote the time size and the time duration the factor m is used, respectively. This treatment is found effective in suppressing the high-frequency noise at the start of the computation. 3.3.2. Fully reflectiÕe boundaries For a fully reflective vertical boundary, the normal velocity must be zero, i.e.: un x q Õn y s 0 or

uPns0

Ž 42 .

in which n is the outward unit normal vector of the reflective boundary and n x and n y are the two components of n along the two spatial coordinates, respectively. The boundary integrals are taken as zero at fully reflective boundaries. In addition to Eq. Ž42., EhrEn s 0 is imposed. 3.3.3. Absorbing boundaries The wave energy arriving at these boundaries from within the fluid domain must be absorbed perfectly. The treatment of this kind of boundary is difficult in numerical modeling. The most commonly used method is the Sommerfeld radiation condition which can be written as: f t q c x f x s 0 and

ft q c y f y s 0

Ž 43 .

where f s u, Õ or h. c x and c y are the phase speeds along the x and y directions, respectively. But in actual computations, especially for irregular waves, it is difficult to know the phase speeds exactly and hence, the boundary can still reflect some wave energy. To eliminate the boundary reflections, a ‘sponge’ layer proposed by Larsen and Dancy Ž1983. is placed in front of an absorbing boundary to absorb the incoming wave energy. On the sponge layer, the surface elevation h and the horizontal velocities Ž u, Õ . are divided by m Ž x, y . after each time step. The factor m Ž x, y . takes the following form after extending the one-dimensional form given by Larsen and Dancy Ž1983. to two dimensions.

mŽ x, y. s

½

exp Ž 2yd rD d y 2yd s rD d . ln a 1

0 F d F ds ds - d

Ž 44 .

in which d is the distance between the point on the ‘sponge’ layer and the boundary, D d is the typical dimension of the elements, d s is the sponge layer width, usually equal to

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

107

one or two wave lengths and a is a constant to be specified. In the following test cases, a s 4 is used and the sponge layer can absorb the arriving waves effectively. The ‘sponge’ layer absorbs nearly all the arriving wave energy. Even if the boundary reflects a small amount of energy, the reflected waves have to pass through the sponge layer again and hence, the boundary behaves as a perfect absorbing boundary. In actual computation, the nodal values of h , u and Õ at the boundary behind the sponge layer can be taken as zero and calculation of the boundary integrals are not required.

4. Test cases 4.1. WaÕe absorption test A perfect absorbing boundary should not allow wave reflections to occur. To test the effectiveness of the sponge layer, several tests are conducted. Fig. 2 shows how the ‘sponge’ layer behaves in a 15-m long wave channel of constant water depth Ž h s 0.5 m. for a solitary wave and a sine wave propagation. The right side of the channel is the absorbing boundary and a 3-m long ‘sponge’ layer is placed in front of it. The computational domain was divided into 200 = 2 elements with D x s D y s 0.075 m and the time step D t was 0.05 s. It can be seen from Fig. 2Ža. that a large portion of the incoming waves can be absorbed. Fig. 2Žb. presents the water surfaces of a sine wave for x ) 9 m at different points of time in one period ŽT s 1.0 s. from t s 49 to 50 s. If reflections had occurred, the reflected waves would have traveled back and forth in the channel about three times. However, the difference between the maximum wave height and the minimum wave height is less than 10%. This means that the effect of wave reflection is slight and the sponge layer is very effective in wave absorption. 4.2. Solitary-waÕe propagation in waÕe channel The simulation of the propagation of a solitary wave over a long distance demonstrates the stability and conservative properties of the numerical scheme. A solitary wave propagation in a wave channel of length 222 m and constant water depth of 0.45 m is simulated. The computation domain was divided into 2960 = 2 elements with D x s D y s 0.075 m and the time step D t was 0.05 s. The incident surface elevation h is:

h s Hsech2

(

(

cs g Ž hqH .

3H 4 h3

Ž x y ct . Ž 45 .

where H is wave amplitude, h is water depth, c is the propagation velocity. Fig. 3 depicts the spatial variations of two solitary waves with amplitudes 0.05 m and 0.15 m at various time steps. The results indicate that the calculated solitary wave profiles are basically identical to the theoretical ones, but the simulated water surface elevations are slightly higher and the difference increases with increasing incident wave height. Also, a

108

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

Fig. 2. Performance of sponge layer.

small dispersive tail is formed behind the evolving solitary wave. These results are similar to those obtained by Wei and Kirby Ž1995. using a finite difference scheme to solve the extended Boussinesq equations derived by Nwogu Ž1993.. Fig. 3 also shows that the numerically predicted phase speed is slightly smaller than the theoretical one, and that the difference increases with increasing wave height. To illustrate this more clearly, Fig. 4 gives the solitary-wave forms at t s 10 s and 60 s which have been translated in position using the theoretical phase speed c. The reason for the reduction of phase speed is that the model equation tends to underestimate the amplitude dispersion.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

109

Fig. 3. Solitary wave evolution in wave channel of constant water depth Ž hs 0.45 m..

In fact, this can be seen by comparing Eq. Ž3. in which b s 1r5 with the solitary wave phase velocity given by Eq. Ž45.. 4.3. WaÕe focusing by a topographical lens Whalin Ž1971. carried out a set of experiments on a wave focusing bottom topography in a wave tank of size 6.10 m = 25.6 m. The topography in the middle portion of the tank was made up of eleven equally spaced semi-circular steps that gradually decreased the water depth from 0.4572 m to 0.1524 m. The following equations describe the topography: 0.4572 0 F x - 10.67 y G h Ž x , y . s 0.4572 q Ž 10.67 y G y x . r25 10.67 y G F x - 18.29 y G Ž 46 . 0.1524 18.29 F x F 21.34 1r2 Ž . w Ž .x where G y s y 6.096 y y , 0 F y F 6.096.

° ¢

~

110

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

Fig. 4. Solitary wave forms at t s10 and 60 s.

Waves with periods T s 1, 2, and 3 s, corresponding to relative depths hrL s 0.306, 0.117 and 0.074 in the h s 0.4572 m portion of the tank, where L is the wave length given by the dispersion relation of the linear theory, were generated and the harmonic amplitudes along the center line of the wave tank were measured at various stations. The

Fig. 5. The bottom topography of computational domain.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

111

Fig. 6. Comparison of calculated and measured wave amplitudes along center line of tank Žsolid lines: the authors; dotted lines: Beji and Nadaoka Ž1996., symbols: experimental data. The lines from top to bottom correspond to the first, second and third harmonics, respectively..

112

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

experimental results have been used by many researchers to verify their numerical models. The bottom topography of the computational domain Ž27.432 = 6.096 m. expressed by Eq. Ž46. is shown in Fig. 5. The domain was divided into 270 = 60 elements for T s 1 s and 180 = 40 elements for T s 2 and 3 s, respectively. At the end of the tank, a 2.432-m long ‘sponge’ layer was used. The time step D t is 0.05 s. The incident wave heights H were 0.039, 0.015 and 0.0136 m for T s 1, 2 and 3 s, respectively. Fig. 6 gives the numerical results, the experimental data and the computed results by Beji and

Fig. 7. Perspective view of the water surface elevation ŽT s 2.0 s, H s 0.015 m and t s 50 s..

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

113

Fig. 8. The computational domain for wave diffraction through a breakwater gap.

Nadaoka Ž1996. using a finite difference scheme. The figure shows that the two sets of computed results agree satisfactorily and are consistent with the experimental results. Fig. 6Ža. is the result for T s 1.0 s. It can be seen from the figure that both the first and second harmonics increase in magnitude in the focal zone. The results are consistent with the experimental results in spite of the scatter in the experimental data. Comparing the two sets of numerical results, the agreement between the first harmonics is slightly better, but the present numerical model underestimates the amplitude of the second harmonic and the oscillations in the wave amplitude that are present in the results of Beji and Nadaoka Ž1996. are absent. The reason may be that the present numerical scheme has a larger damping effect than that of Beji and Nadaoka Ž1996. due to the method used to calculate the derivatives as given by Eq. Ž30.. This aspect will be studied further. For T s 2.0 s, Fig. 6Žb. shows that the model gives a good prediction of the first, second and third harmonic wave amplitudes. The two sets of numerical results are close to each other. The oscillations in the computed wave amplitude, which are now of much

Fig. 9. Sketch of the finite elements near the breakwater tip.

114

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

longer wave length, are also present, but the wave lengths are slightly different from that of Beji and Nadaoka Ž1996.. This is also probably due to the damping introduced in the present scheme as mentioned above. Fig. 6Žc. for T s 3.0 s shows that the numerical model overestimates the amplitude of the first harmonic and underestimates the second and third harmonic amplitudes. The

Fig. 10. Comparison of the computed and experimental diffraction coefficient contours Žsolid line: numerical, dashed line: experimental..

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

115

Fig. 10. Žcontinued..

results are also similar to that obtained by Liu et al. Ž1985. who explained that the differences between the numerical and experimental results were probably due to the effect of viscous damping and the relatively short evolution distance. Another possibility may be due to the presence of reflected free waves, leading to energy exchanges

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

116

between harmonics. In general, the agreement between the two sets of numerical results is better than that for T s 1 and 2 s. A perspective view of the fully-developed water surface elevation at t s 50 s is shown in Fig. 7 to illustrate the wave focusing phenomenon. 4.4. WaÕe diffraction through a breakwater gap To further evaluate the applicability of the numerical model in real situations, wave diffraction through a breakwater gap is simulated and the results are compared with laboratory data. Laboratory experiments on wave transformation through a breakwater gap were carried out by Yu et al. Ž1998.. In their experiments, a breakwater of thickness 0.35 m was parallel to and located 7 m in front of the wavemaker. The gap in the center was formed by two semi-circular tips. The experimental area was 26 m = 27 m. Two gap widths, B s 3.92 m and 7.85 m, were studied. Wave absorbers of thickness 0.8 m were placed on all the three sides of the experimental area and in front of the breakwater except in the vicinity of the breakwater tips to eliminate wave reflections. Regular, uni-directional and multi-directional wave experiments with principal wave directions of 908, 758, 608 and 458 were conducted. The wave height or significant wave height was 0.05 m and the wave period or peak wave period was either 0.92 s or 1.20 s. The water depth was constant at 0.4 m for wave diffraction. More detailed information about the laboratory experiments is given in Yu et al. Ž1998.. From the measurements, the diffraction coefficient K d is computed using the following equation: Kd s

Hd Hi

Ž 47 .

where Hd is the diffraction wave height for regular waves or significant wave height for irregular waves and Hi is the incident wave height, respectively. In this paper, the diffraction of regular waves through a breakwater gap corresponding to the experiments is modeled. Because of the long computation time, the computed domain is only 24 m = 16.4 m as shown in Fig. 8 which is smaller than the experimental area. The incident boundary is 4 m in front of the breakwater. For B s 3.92 m and 7.85 m, 17985 and 18069 four-node elements, corresponding to 17584 and 17696 nodes, respectively, are used. The semi-circular breakwater tip is approximated by four line segments as shown in Fig. 9. The time step D t s 0.05 s and the total time of computation is 40 s. After the calculation, the computed wave diffraction coefficients are computed using Eq. Ž47.. Fig. 10 is the computed wave diffraction diagram together with the experimental results for directions u 0 s 908 and 458. Fig. 11 gives the cross-section distribution of the diffraction coefficients at YrL s 3. The figures show that the values of the computed diffraction coefficient are basically identical to the experimental values, but the experimental results are in general a little higher than those computed. This may be due to the difference between the absorbing characteristics of the sponge layers used in the numerical model and those of the physical absorbers, at the absorbing boundaries, especially those located in front of the breakwater. The width of the sponge layer is

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

Fig. 11. Cross-sectional distribution of the diffraction coefficient at Y r Ls 3.

117

118

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

Fig. 12. Perspectives of water surface at t s 40 s.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

119

about one wavelength and is thicker than the thickness of the physical absorbers Žs 0.8 m. used in the experiments. The sponge layer will affect the wave propagation behind the breakwater. In addition, the semi-circular breakwater tips are approximated by four line segments. The slight difference in the shape of the breakwater tip used in the numerical calculation has an effect on the reflection of waves from the breakwater tips. Fig. 12 gives the perspective of the wave surface at the instant t s 40 s. The figure also shows the effect of reflection by the breakwater tips on the waves before and after the breakwater. Especially for oblique incident waves, the waves in front of the breakwater are apparently affected by reflected waves because the incident boundary does not allow the absorption of reflected waves from the computational domain. This aspect should be studied further.

Fig. 13. Computational domain and finite element grid.

120

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

Fig. 14. Comparison of computed, experimental and linear theory runup profiles.

4.5. WaÕe runup around circular cylinder Finally, the advantage of the finite element method is demonstrated by simulating the wave runup around a circular cylinder for which the finite difference method is difficult to be applied. The physical problem is illustrated in Fig. 13Ža.. The length and width of the computational domain are 6.4 m and 3.4 m, respectively. A total of 4964 elements were used as shown in Fig. 13Žb.. The water depth h is 0.08 m and the radius of the cylinder a is 0.25 m. The incident wave height Hi and wave number k are 0.03 m and 5.2 my1 , respectively resulting in ka s 1.3 and l s Hirk 2 h 3 s 2.14. The parameter l is equivalent to the Ursell number Žs 4p 2rl2 .. The incident wave period T is 1.4 s and the time step D t s 0.04 s. Isaacson Ž1978. performed laboratory experiments on wave runup around a circular cylinder in the shallow water range. The values of ka and l chosen for the present simulation are the same as those in one of the experiments carried out by Isaacson Ž1978.. Fig. 14 is the computed wave runup profile RŽ u .rHi around the cylinder together with the experimental results as well as theoretical predictions from linear wave theory. It can be seen that the linear theory underestimates the runup and the numerical results agree better with the experimental results, especially the maximum runup at u s 1808.

Fig. 15. Perspective of water surface at t s10 s.

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

121

The differences between the numerical and experimental results are mainly due to reflected waves from the cylinder being reflected again by the incident wave boundary in the computation. The contamination of the numerical results by small reflected waves can be seen from Fig. 15 which gives the perspective of the wave surface at the instant t s 10 s. The wave crest varies slightly in the y direction near the incident wave boundary Ž x s 0.. As mentioned in previous sections, the problem of wave reflection by the incident boundary requires further study.

5. Conclusions A numerical model based on the improved Boussinesq equations using the finite element method is developed in this paper. The time integration is performed using the fourth-order Adams–Bashforth–Moulton predictor and corrector scheme. Linear quadrilateral elements are used and the first spatial derivatives of the water surface elevation at any node are determined by the weighted average of the corresponding derivatives in elements surrounding the node. This corresponds to the central difference approximation of the first derivatives in case of rectangular elements. The outgoing waves are effectively absorbed using a ‘sponge’ layer which can be easily implemented in the model. To save computer storage, only the non-zero elements are stored and the global equations are solved by iteration. The numerical results for solitary wave propagation over a long distance demonstrate that the numerical model is stable and conservative. The topographical lens problem is modeled and the results show that the model can predict the expected harmonics generated by the lens topography and are in agreement with the experimental results given by Whalin Ž1971. in spite of some relatively minor discrepancies. For wave diffraction through a breakwater gap, the model provides satisfactory diffraction diagrams and the results are consistent with the experimental data of Yu et al. Ž1998.. The slight difference is basically due to the fact that the numerical model does not correspond exactly to the physical simulation because of limitations of the model. The final test case of wave runup around a circular cylinder further demonstrates the versatility of the proposed model. In general, it can be concluded that the present model is capable of providing satisfactory results in engineering applications. The diffraction of waves through a thick breakwater gap with semi-circular tips and the wave runup around a circular cylinder are two such examples.

Acknowledgements The work reported in this paper is supported by both the Hong Kong Research Grants Council ŽNo. PolyU 57r96E. and the National Natural Science Foundation of China ŽContract No. 59479006..

122

Y.S. Li et al.r Coastal Engineering 37 (1999) 97–122

References Abbott, M.B., McCowan, A.D., Warren, I.R., 1984. Accuracy of short wave numerical model. J. Hydr. Eng. 110 Ž10., 1287–1301. Antunes do Carmo, J.S., Seabra Santos, F.J., 1996. On breaking waves and wave–current interaction in shallow water: a 2DH finite element model. Int. J. Numer. Methods Fluids 22, 429–444. Antunes do Carmo, J.S., Seabra Santos, F.J., Barthelemy, E., 1993. Surface waves propagation in shallow ´ water: a finite element model. Int. J. Numer. Methods Fluids 16, 447–459. Beji, S., Nadaoka, K., 1996. A formal derivation and numerical modeling of the improved Boussinesq equations for varying depth. Ocean Eng. 23 Ž8., 691–704. Dingemans, M., 1997. Water wave propagation over uneven bottoms. Advanced Series on Ocean Engineering, Vol. 13, Parts 1 and 2. World Scientific, Singapore. Hamm, L., Madsen, P.A., Peregrine, D.H., 1993. Wave transformation in the nearshore zone: a review. Coastal Eng. 21, 5–39. Howard, D., Connolley, W.M., Rollett, J.S., 1990. Unsymmetric conjugate gradient methods and sparse direct methods in finite element flow simulation. Int. J. Numer. Methods Fluids 10, 925–945. Hughes, T.J.R., 1987. The Finite Element Method. Prentice-Hall, Englewood Cliffs, NJ. Isaacson, M. de St.Q., 1978. Wave runup around large circular cylinder. J. Wtrwy., Port, Coast. Oc. Division, ASCE 104 Ž1., 69–79. Kashiyama, K., Ito, H., 1995. Three-step explicit finite element computation of shallow water flows on a massively parallel computer. Int. J. Numer. Methods Fluids 21, 885–900. Kirby, J.T., 1990. Modeling Shoaling Directional Spectra in Shallow Water. Proc. 22nd Int. Conf. Coast. Eng. ASCE, New York, NY, pp. 109–121. Larsen, J., Dancy, H., 1983. Open boundaries in short wave simulation—a new approach. Coastal Eng. 7 Ž3., 285–297. Li, Y.S., Zhan, J.M., 1998. A 3-dimensional finite element method for stratified coastal seas. J. Hydr. Eng. ASCE 124 Ž7., 699–703. Liu, P.L.-F., Yoon, S.B., Kirby, J.T., 1985. Nonlinear refraction–diffraction of waves in shallow water. J. Fluid Mech. 153, 185–201. Madsen, P.A., Sørensen, O.R., 1992. A new form of the Boussinesq equations with improved linear dispersion characteristics: Part 2. A slowly-varying bathymetry. Coastal Eng. 18, 183–204. Madsen, P.A., Sørensen, O.R., 1993. Bound waves and triad interactions in shallow water. Ocean Eng. 20 Ž4., 359–388. Madsen, P.A., Murray, R., Sørensen, O.R., 1991. A new form of the Boussinesq equations with improved linear dispersion characteristics. Coastal Eng. 15 Ž4., 371–388. Mase, H., Kirby, J.T., 1992. Hybrid Frequency-domain KdV Equation for Random Wave Transformation. Proc. 23rd Int. Conf. Coast. Eng., ASCE, Venice, pp. 474–487. Nwogu, O., 1993. An alternative form of the Boussinesq equations for nearshore wave propagation. J. Wtrwy, Port, Coast. Oc. Eng. ASCE 119 Ž6., 618–638. Peregrine, D.H., 1967. Long waves on a beach. J. Fluid Mech. 27, 815–827. Schaffer, H.A., Madsen, P.A., 1998. Discussion of ‘a formal derivation and numerical modelling of the ¨ improved Boussinesq equations for varying depth’. Ocean Eng. 25 Ž6., 497–500. Wei, G., Kirby, T., 1995. Time-dependent numerical code for extended Boussinesq equations. J. Wtrwy, Port, Coast. Oc. Eng. ASCE 121 Ž5., 251–261. Whalin, R.W., 1971. The limit of applicability of linear wave refraction theory in a convergence zone. Res. Rep. H-71-3, U.S. Army Corps of Engrs., Waterways Expt. Station, Vicksburg, MS. Witting, J.M., 1984. A unified model for the evolution of nonlinear water waves. J. Comput. Phys. 56, 203–236. Yu, Y.X., Liu, S.X., Li, Y.S., Wai, W.H., 1998. Refraction and diffraction of random waves through breakwater. Submitted to Ocean Engineering.