The modelling of infinitesimal amplitude water waves using an improved boundary element method

The modelling of infinitesimal amplitude water waves using an improved boundary element method

002&7225/89 $3.00+ 0.00 Inr. 1. Engng Sci. Vol. 27, No. 10, pp. 1217-1222,1989 Printed in Great Britain. All rights reserved Copyright @ 1989Pergamo...

404KB Sizes 2 Downloads 70 Views

002&7225/89 $3.00+ 0.00

Inr. 1. Engng Sci. Vol. 27, No. 10, pp. 1217-1222,1989 Printed in Great Britain. All rights reserved

Copyright @ 1989Pergamon Press plc

THE MODELLING OF INFINITESIMAL AMPLITUDE WATER WAVES USING AN IMPROVED BOUNDARY ELEMENT METHOD T. T. AL-SHEMMERI’

and 0. M. S. I-LAMED2

‘Department ‘Department

of Mechanical and Computer Aided Engineering, Staffordshire Polytechnic, Beaconside, Stafford ST18 OAD, U.K. of Mathematics/College of Science, King Saud University, Riyadh, Saudi Arabia

Abstract-A description of the improved boundary element method is presented. Its application to predict the motion of free surface flows is described, and a model of an infinitesimal water waves is demonstrated in particular.

1. INTRODUCTION

The modelling of the dynamic behaviour of free surface flows is one of the most challenging problems in computational fluid dynamics (CFD). Difficulties are accentuated by the location of the free surface which is not known a priori. The non-linearity and transient nature of the system is an added complication to the problem. Only the most general solution techniques which impose no restrictions on the shape of the boundaries offer the flexibility required to model such a problem. Furthermore, the fine discretization needed to resolve highly concentrated potential fields results in large computer expenditures which must be minimized when and wherever possible. In these situations modelling techniques such as finite difference or finite element methods, become less efficient as it is required to discretize the whole domain. In contrast, the boundary element method (BEM) which requires only to discretize the boundary offers a drastic reduction in the computations required. Furthermore, there is no real need to solve at internal points within a free surface flow as the name implies, one is only interested in the movement of the free surface itself. However, the conventional BEM has its pitfalls. Errors due to singularities become unacceptable, particularly with the transient problem, as these errors propagate and eventually cause divergence. This study shows that improvements in the general techniques can overcome such problems and provide accurate solutions. The method of overcoming such errors does not contribute to an increase in CPU time and memory since this method of curing singularities is achieved by an innovative approximation of the integrals at positions of discontinuity. The purpose of this paper is to give a brief description of the modifications leading from a conventional BEM approach to a potential flow to the procedure used in this study, and to present the application of the modified approach to the modelling of water waves and associated problems. Full analysis of a given test case is presented.

2.

FORMULATION

OF

FREE

SURFACE

ELEMENT

FLOW

USING

BOUNDARY

METHOD

Free surface flows are governed by the energy equation of a non-steady given by

incompressible

fluid

Coupled with

(2) 1217

1218

T. T. AL-SHEMMERI

and 0.

M. S. HAMED

which represents the instantaneous free surface elevation f in terms of the velocity potential @. Rewriting equations (1) and (2) into finite time difference form using weighting factors 6, and e2 to specify solution scheme. In the present paper the Crank-Nicholson scheme was chosen (i.e. 0, = t& = 0.5). Hence @K+l

-aK

[ 1 IVQT~~K + e+ =2

At

+F)“”

+ (I-

e,)(gc +:)K]

(3)

and c

Kc1 -

At

5;”

= [ e,($)K+l

+ (1 - e@JK]/cOs

PK.

Superscripts K and K + 1 belong to time steps t and t + At, respectively. between equations (3) and (4) yield Q,k+l =

a& -g(~t)ze,e~ +

[a@/anlk+’ [cos/3]” -

(4) Elimination

of cK”

At{: Kp)z +

(s321K

gck+ At(l - ~,)~,g[wan]V[c0s 01”+ h[p,iplk+l + (1 - e,)[p,/p]“)j.

(5)

It is well known [l] that the potential @ satisfies the following integral equation

where @*=&ln

1 ,

(7)

0r is

the fundamental solution of a 2-D In order to solve equation (6) for boundary. The procedure is thus to discretize and d@/dn is assumed to be linear. surface Q, @ or a@/& is represented

isotropic case. ai, the values of @ or a@,/& are required on the entire the surface 52 into boundary elements. The variation of Q, That is, over each interval (boundary element) along the by

J = 1 to N, @=(l-

PP(4j)

+PWqj+d I

da/an = (1 - p) E(e)+P$!(qj+l) an

09

I’

where qj and qj+i are the end points of dQ, and ,Uis a linear function non-dimensionalized such that its value lies between 0 at qj and 1 at qj+l+ By choosing nodal points in succession and numerically integrating equation (6), then applying the appropriate boundary conditions, a system of N x N equations is formed such that

Rearranging (9) to set all known arrays in the RHS vector and all unknowns on the LHS vector we thus obtain a system matrix of the following form Ax=b,

which is solved by the Gaussian elimination procedure.

WV

Once the solution

is obtained

on the free surface, QK+l can be calculated using equation (5) and also the new surface elevation p”+’ using equation (4). The solution then proceeds to the next time step and so on. REMARK 1. A well-known problem arises at corners where the true value of da/an undergoes a sudden change caused by the change in direction of the normal at the boundary. The classic

Modelling water waves using an improved boundary element method

1219

way of dealing with this is to replace the single node by two nodes a short distance apart [2]. However, this only provides a limited improvement. A very elaborate method [l] was introduced by the author which consists of allowing the two nodes to move together and coalesce, and then by recognizing that code only elements adjacent to a corner are integrated analytically. This provided a dramatic improvement in accuracy. A test with infinitesimal wave theory proved this treatment of corners, with errors of the order of errors in numerical integration of far field elements. REMARK2. The second feature of this paper is the adaptive integration procedure based on the Gauss-Legendre polynomials. A q-point Gauss-Legendre quadrature rule is used to evaluate the approximate integrals. The error incurred when such a rule is used is

Eq =[(2q

(404 + 1)((2q)!)3]

h24+1f24(P)9

for some p E (q, b),

see Stroud and Secrest [3]. Now a/an log(l/n) = cos 8/R and the 2qth derivative of the integrand can be estimated from dz9 cos 0 *

( y-

(2q)! ++I >

R being the distance of the fixed point P(x, y) from the centre point & (x, y) on the integration

interval. The error E, can therefore be estimated u priori and q chosen so that E4 < hj s/2 where E is a measure of the accuracy required in the final solution and hj is the jth element length. Details of the error analysis are given in [4]. 2.1 Standing waves of infinitesimal amplitude Analytically, only when

the simplest case of water waves is the linear wave theory. This theory is valid (a/A)[A/h]”

is small, i.e. for waves of small amplitude.

A standing wave (Fig. 1) travelling from left to right in a tank can be defined in terms of its spatial location (x) and time (t) as 5 = a sin(kx

(11)

where 5 is the local free surface height at distance x from the source of the wave origin. k is the wave number and (J is the reciprocal of the wave length. A standing wave of the form given in equation (11) can be generated by three possible methods: (a) assuming a starting elevation p on free surface;

/ Y / 4 / If / 47

/ /) / / / N,’ / / 1 ““‘r”“““““““““” H

L

Fig. 1. Definition of free surface geometry.

*

h

1220

T. T. AL-SHEMMERI

and 0.

M. S. HAMED

(b) Assuming known velocity on the free surface by differentiating

aa

n

(c) assuming a pressure distribution

(11)

(12)

= aa sin(kx)cos(of); on the free surface given in [5] as P = -y

sin(b).

(13)

Therefore, three alternative test cases are available to check the applicability approach to model water waves.

of the present

2.2 Description of the test problem The problem of particular interest seen in Fig. 1 consists of a tank 1 m long and 0.5 m deep. The amplitude was chosen as 1 cm for which the Ursell number should be less than unity (14) in order to ensure that infinitesimal amplitude linear wave theory is valid. The circular frequency for this wave is u = qgk tank(kh) = 5.3166, and the time period is TP= (2~r/a) = 1.1818 s. The bottom comers suffer no discontinuity and therefore only free surface corners were treated for singularity [l]. The total number of nodes used are 27, 12 of which are on the free surface. In solving the problem, after each time step the elements on the sides of tank were updated in accordance with the current location of the free surface. In particular, the end nodes (i.e. the crest and trough) were adjusted so that the nodes on the sides were maintained with equidistant elements. 2.3 Quadratic Interpolation of Wave-Crests and Troughs The computed results obtained by a numerical method do not necessarily produce values of elevation at the exact locations of crests or troughs. This is due to the finite time step and geometrical discretization. Therefore it is necessary to locate the wave-crests and troughs approximately using the numerical computed results. The approximate times of the crests and troughs are readily identified from the numerical data and three points were chosen to span each one. These points were fitted to a quadratic equation.

5‘= a + bt + ct’,

(19

where a, b and c are constants. For constant time intervals, it is readily shown that the maximum (or minimum elevation) given by g,,, =

c2

is

(L - f;l)’

WC3

-

X2

+

5‘1)’

at time

(53- L-W f=t2-2(g3-2C2+

cl)’

(17)

3. RESULTS The three initial conditions described earlier were applied to obtain solutions for the standing wave problem of infinitesimal amplitude (0.01 m amplitude in a tank of 0.5 m depth and wavelength of 2 m) with N = 25. Results appear in Table 1 and demonstrate the inaccuracy of solutions without singularity treatment.

Modelling water waves using an improved boundary element method Table 1. BEM results of standing waves. Linear-untreated (I = 0.01 m

Method Elevation

1221

case N = 25, 7’S = 0.02 s,

Crest’s height

Relative error at crest

Trough’s height

Relative error at trough

0.011906 0.010815 0.011641 0.010522

0.1906 0.0815 0.1641 0.0522

0.011716 0.010787 0.011785 0.010032

-0.1716 -0.0787 -0.1785 -0.0132

av 0.1221

av -0.1105

Relative error in wave period after complete period is +lO.%% Velocity

0.010776 0.011111 0.008360 0.011842

0.0776 0.1111 -0.1640 -0.1842 av

0.010603 0.010774 0.088656 0.011588

0.052225

-0.0603 -0.0774 +0.1344 -0.1588 av -0.040525

Relative error 7.84% in period Pressure

0.010832 0.011022 0.008382 0.011477

0.0832 0.1022 -0.1618 0.1477 av

0.010643 0.010695 0.008688 0.011721

0.042825

-0.0643 -0.0695 +0.1312 -0.1721 av -0.043675

Relative error = 9.79% in period

The error in the computed wave period was 10.96% for the elevation start, 7.84% for the velocity start, and 9.79% for the pressure start. The average error of computed crests was 12.21, 5.22 and 4.28%, respectively and of troughs was 11.05, 4.05 and 4.36%) respectively. Table 1 indicates that the velocity start, with an initial horizontal free surface, appeared to be relatively the best of the three approaches for this particular test case. Table 2 contains results obtained with singularity at corners being correctly treated. A glance at these results immediately demonstrates the drastic improvements in these solutions. Errors in elevation of the wave at the crests were only 0.0725,0.4425 and 0.1675% for the three initial conditions of elevation start, velocity start and pressure start, respectively and 0.41, 0.07 and 0.3425% at the troughs. Errors in wave periods computed were 0.335, 0.363 and 0.448%, respectively. Thus the results were dramatically improved and became relatively accurate compared with linear wave theory. Table 2. BEM results of standina waves: Treated case. N = 27. T.9 = 0.02 s. a= 0.01 m

Method Elevation

Crest’s height 0.01081 O.Olw22 0.089999 O.W9998

Relative error at crest 0.001 0.0822 -O.oool -0.wJ2

Trough’s height

Relative error at trough

0.009964 0.009927 0.009960 0.009985

0.0036 0.0073 0.0040 0.0015

av 0.00725

Velocity

Relative error in wave period 0.335% 0.010015 0.0015 0.089973 0.010072 0.0872 0.009997 0.010880 0.0080 0.010021 0.010010 0.0010 0.009981 av 0.004425

Pressure

av 0.0041

Relative error in wave period 0.363% 0.00999 -0.001 0.009948 0.010041 0.0041 0.009973 0.010050 0.0050 0.009991 0.@39986 -0.0014 0.009951 av 0.001675 Relative error in wave period = 0.448%

av0.0007

0.0052 0.0027 O.OCW 0.0049 av 0.003425

1222

T. T. AL-SHEMMERI 4.

and 0. M. S. HAMED

CONCLUSIONS

The accuracy stability of the BEM approach to time-varying water demonstrated by comparisons with known solutions for standing waves.

waves has been

REFERENCES [l] T. T. AL-SHEMMERI, Ph.D. Thesis, UMIST, U.K. (1983). [2] C. A. BREBBIA, The Boundary Element Method for Engineers. Pentech Press, London (1978). [3] $9E),STROUD and D. SECREST, Gaussian Quadrature Formulas. Prentice Hall, Englewood Cliffs, New Jersey [4] 0. M. S. HAMED, Ph.D. Thesis, University of Salford, U.K. (1982). [S] H. LAMB, Hydrodynamics, 6th Edn. Cambridge University Press (1932). (Received 4 January 1989)

NOMENCLATURE a

; k K

P,

t

T

x

wave amplitude acceleration due to gravity water depth wave number time step iteration counter pressure term time wave period horizontal coordinates

wavelength angle of free surface with horizontal boundary surface surface elevation weighting factors density wave property potential function domain