Эффективный численный метод моделирования сильно нелинейных волн Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 3, № 3, 1998

AN EFFICIENT NUMERICAL METHOD FOR CALCULATION OF STRONGLY NONLINEAR

WATER WAVES

B. E. PROTOPOPOV Lavrentyev Institute of Hydrodynamics, Novosibirsk, Russia e-mail: [email protected]

Предлагается эффективный численный метод моделирования сильно нелинейных волн, включая процесс обрушения. Метод основывается на смешанном эйлерово— лагранжевом подходе. Основная идея этого подхода распространяется на всю область течения. Это означает, что не только свободная поверхность, но и жесткие границы и внутренние узлы сетки движутся как частицы жидкости. Это позволяет лучше визуализировать течение, генерировать сетку только одни раз и влечет экономию памяти и времени счета. Еще одно отличительное свойство предлагаемого численного метода заключается в том, что дискретизация по пространству проводится с помощью сплайнов, при этом производные функций вычисляются аналитически, дифференцированием соответствующий аппроксимирующих сплайнов. Численный метод протестирован на трех задачах: эллипс жидкости***, обрушивающаяся *** и возмущенная уединенная волна большой амплитуды. Результаты тестовых вычислений показывают высокую точность и широкую применимость предлагаемого метода.

1. Introduction

The aim of this work is to develop an efficient numerical method which is capable of treating strongly nonlinear water waves, including overturning ones.

When someone intends to simulate overturning waves, he has, first of all, to formulate the problem with parametric representation of a free water surface. Such formulation looks, in general, more complicated, compared to the "traditional" one (when free-surface parameters coincide with the horizontal Eulerian coordinates). But the complexities can be effectively removed by choosing free-surface parameters to be conserved in fluid particles or, equivalently, to coincide with the Lagrangian coordinates. This idea, known as the mixed Eulerian — Lagran-gian (MEL) formulation, was first successfully realized by Longuet — Higgins & Cokelet [1] and then intensively exploited by many investigators (let us name only a few — [2-12]).

Since the MEL formulation implies the mobility of free-surface nodes, a computational grid has to be generated at each time step. As a result, the requirements of CPU time and memory for grid generation become significant. To diminish these requirements, most of the investigators [1-7] has applied the boundary-integral equation (BIE) method: the problem is reformulated on the fluid-domain boundary so that the dimension of a computational grid is

© B.E. Protopopov, 1998.

reduced by unity. This advantage of the BIE method is, however, essentially suppressed by its shortcoming: matrices of the resultant linear-equation systems are always dense. This makes it reasonable to develop alternative, space-discretization (SD) methods, which though require a computational grid in the whole fluid domain but usually yield banded matrices. Examples of successful application of such methods, namely the finite-difference and finite-element ones, can be found in [8-10] and [11, 12], respectively.

The basic distinctive feature of the present numerical method, which belongs to the SD ones, is that a computational grid is generated only once, at the initial moment of time. Then all grid nodes are moved as fluid particles. Such approach saves the CPU time and memory and makes water motion more visual. But care should be taken to keep a computational grid smooth. Here this care is provided mainly by application of splines for space discretization (this is one more distinctive feature of the method). Derivatives of any function with respect to space coordinates are evaluated analytically, by differentiating the spline which approximates this function. The splines are chosen to be of second order as seeming most optimal, i.e. providing high enough accuracy with low expense. The present numerical method possesses some other peculiarities which, however, are of secondary significance and can be seen below, in the method description.

The general formulation of the problem on unsteady motion of incompressible inviscid fluid with a free surface is given in the next paragraph in terms of the generalized (Lagrangian) coordinates. Paragraph 3 is devoted to the detailed description of the numerical method developed. Paragraph 4 contains the results of test calculations with some discussion of these results.

2. Governing equations

Only the two-dimensional problem is considered here. Being accepted for simplicity, this limitation is, however, non-principal: the method proposed can be easily extended to the three-dimensional case.

Let x be a horizontal coordinate in the Cartesian system, z a vertical one, and t time. The fluid is bounded by a free surface and by rigid boundaries, including a bottom, left and right walls. All the rigid boundaries are prescribed in a general form:

the bottom : B(x,z,t) = 0,

the left wall : L(x, z, t) = 0,

the right wall : R(x, z, t) = 0, (1)

with the functions B(x,z,t), L(x,z,t), R(x,z,t) being differentiable with respect to all their arguments.

2.1. The interior equations

The fluid is assumed to be incompressible, inviscid and moving irrotationally. Under the assumptions accepted, the velocity vector possesses a potential $(x,z,t) which is a harmonic function inside the domain (hereinafter letter subscripts indicate partial differentiation):

Before writing appropriate boundary conditions, we first introduce new (generalized) space coordinates, £ and (, in order to match the fluid domain onto a fixed rectangle, say [£0, £M]x x [Z0, (N], and thus to remove some difficulties associated with the implementation of boundary conditions. In terms of the generalized coordinates, Laplace's equation becomes a general elliptic equation which is written here in the divergent form

uç + wz = 0, (2)

with

u = (Y0i + p0z )M w = (a0z + P0ç )M (3)

J = xç zz — zç xz, (4)

a = x| + z|, P = — (xç xz + zç zz), y = x^ + z^. (5)

All the quantities to be determined, namely x, z, 0, u, w, J, a, P, y, are considered now as functions of Z and t. In particular,

0(£,Z,t) = $(x(£,Z,t),z(£,Z,t),t).

From the system (2)-(5), we have the following restrictions on the coordinate transformation: all the derivatives xç, xz, zç, zz must exist and the Jacobian J must be different from zero. The problem of constructing the coordinate transformation with the properties just mentioned, as well as some other desirable properties, is referred to as the grid generation problem. Here a computational grid is generated only once, at the initial instant of time:

x = xo(£,Z ), z = zo(£,Z ) (t = 0). (6)

Usually the fluid domain is not complicated initially, so that some simple algebraic procedure can be applied to generate a computational grid. Then all grid nodes are moved as fluid particles (in other words, the generalized coordinates coincide with the Lagrangian ones):

xt = zt = , (7)

or, in terms of u and w:

xt = (xç U + xz w)/J, zt = (zç U + zz w)/J. (8)

Resolving the system (8) with respect to u and w, one can obtain

U = zz xt — xz zt, w = — zç xt + xç zt. (9)

From (4) and (9), it is easy to verify the quite obvious fact: the fluid incompressibility implies the conservation of elementary fluid volumes, that is the Jacobian invariance:

Jt = uç + wz = 0.

This is one of the most attractive benefits of use of the Lagrangian coordinates.

2.2. The boundary equations

Since the conditions on all the rigid boundaries are quite similar, we consider only the bottom which is matched onto the line Z = Z0 in the (£, Z)-plane. The bottom equation (1) is equivalent to the following partial differential equation

Bx xt + Bz zt = - Bt, (10)

with initial data

B(xo,zo, 0) = 0. (11)

Differentiation of the equation (1) along the bottom, that is with respect to £, yields

Bx x, + Bz z, = 0. (12)

The rigid-boundary conditions are kinematic ones, expressing that at any moment of time the domain boundaries consist of the same fluid particles: if fluid-particle coordinates satisfy (11) initially, then fluid-particle velocities must satisfy (10). By substituting (8) into (10) and then using (12), we obtain the following boundary (bottom) condition for the elliptic system

(2)-(5):

w = - B x + B z Bt (Z = Zo). (13)

Bx Xz + Bz Zz

Now the relations (8), with w given by (13), do not guarantee the motion of fluid particles strictly along the bottom, since in deriving (13) we have used the condition (12) which is a consequence from (1), i.e. from (10), (11). So the equations governing the motion of fluid particles along the bottom should be deduced directly from (10) and the first of (9):

-Bt xz + Bz u Bt Zz + Bx u

Xt = R x + R z ' Zt = - R x + R z (Z = Z0)- (14)

Bx xz + Bz Zz Bx xz + Bz Zz

The boundary values of u and the equations of motion of the side-wall nodes can be obtained in the quite similar manner:

u = - OT^ L = ^ (15)

xt = - Lf+rf■ Zt = t'Ztt* « = £o); (")

u = - Rxx, + r,z, R « = )- (17>

x = - Rt x, + R, w z = -Rt z, + Rx w = m) (18)

Rx x, + Rz z,' Rx x, + Rz z,

The free-surface condition is dynamic one, expressing the pressure continuity across a free surface. This is Bernoulli's equation written in terms of the Lagrangian coordinates (g is the gravitational acceleration):

& = 2(x? + Zt2) - gz (Z = Zn), (19)

where the expressions for xt and zt are the same as for inner nodes:

xt = (x, u + xz w)/J, zt = (z, u + Zz w)/J (Z = ZN)• (20)

The equation (19) is of evolutionary type and, hence, needs appropriate initial data:

& = &o(£) (Z = Zn, t = 0). (21)

2.3. The corner equations

In the left-wall-bottom corner, except (10), (12), the similar conditions, but with the left-wall function L(x, z,t), are valid:

Lx xt + Lz zt = — Lt, (22)

Lx xz + Lz zz = 0. (23)

The relations (10) and (22) immediately give us xt and zt. Substitution of these derivatives into (9) leads to the expressions for u and w, which can be then simplified using (12) and (23). Thus, the corner values of velocities and the equations of motion of the corner nodes look as

u = - fff + Bz L L (i = €□, C = Co), (24)

Lx Bz Bx Lz

w = - LLx?+LzL Bt (C = Co, C = Co), (25)

* ^x ; z ; x * ^z

Lz Bt — B7 Lt Bx Lt — Lx Bt /x . ,

xt = L z B^ B L ' zt = L xB^ B L (C = Co, C = Co). (26)

Lx Bz Bx Lz Lx Bz Bx Lz

Similarly for the right-wall-bottom corner:

U = - Bff + Bz Rz Rt (C = Cm, C = Co), (27)

x z x z

W = - Rff + RRz !! Bt (C = Cm, C = Co), (28)

x z x z

Rz Bt - Bz Rt Bx Rt - Rx Bt

xt = R B _ B R , zt = R B _ B R (C = CM, C = Co). (29)

x z x z x z x z

The equations which are valid in the left-wall-surface corner are combined from (15), (19) and (16) to yield:

u = - L x + L z Lt (C = Co, C = Cn), (30)

Lx x5 + Lz z5

0t = 1 (xt2 + zt2) - gz (C = Co, C = Cn), (31)

Lt x? + Lz w -Lt z5 + Lx w t . . ,

xt = T-—7-, zt = --—- (C = ^ C = CN). (32)

Lx x? + Lz z? Lx x? + Lz z?

Similarly for the right-wall-surface corner:

U = - R-Rt (C = Cm, C = Cn), (33)

Rx x? + Rz z?

0t = 1(x2 + z2) - gz (C = Cm, C = Cn), (34)

xt = - Rtxg + Rzw, zt = -Rtzg + Rxw (C = CM, C = CN). (35)

Rx x^ + Rz z^ Rx x^ + Rz z^ '

3. Numerical method

In the previous paragraph, the problem under investigation is formulated as the system of partial differential (first-order) equations (2)-(6), (8), (13)-(21), (24)-(35). Being very complicated, this system is solved numerically. The following three basic parts of the numerical method can be distinguished: (i) discretization in time, (ii) discretization in space, (iii) solution of the elliptic problem.

3.1. The time-stepping procedure

To demonstrate the time-stepping procedure developed, let us consider the model equation

Pt = q, (36)

where q does not depend on time-derivatives of p. Discretization in time is carried out according to the Crank—Nicolson (second-order) scheme:

pn+1 _ pn 1

P_P _ 1 ( n+1 , n\

T 2

or, equivalently:

(q 1 + qn),

pn+1 - 2 qn+1 = rn+1. (37)

Here t is a time step, the superscripts indicate time-step numbers, the right-hand side of (37) contains all the quantities known from the n-th time level:

rn+1 = pn + 2 qn

Using (37) written with the superscript n, one can replace the last relation by

rn+1 = - rn + 2 pn, (38)

which gives a simple rule for determining the right-hand side of equation (37), in a successive manner, starting from

r0 = p0 - 2 q0, (39)

with p0 and q0 being known from initial data.

Since generally the dependence of q on p is nonlinear, the equation (37) is calculated itera-tively (the second superscripts indicate iteration stages):

n+1,k+1 _ n+1 | T n+1

p r + q

starting from

pn+1,fc+1 = rn+1 + _ qn+1,fc (40)

pn+1,0 = rn+1 + 2 qn or, equivalently,

pn+1>0 = - pn + 2 rn+1. (41)

These four rather simple equations, namely (38)-(41), comprise the time-stepping procedure. This procedure is coupled with iteration of the elliptic problem, and the interruption criteria is common for both the procedures (see below).

All the evolutionary equations from the previous paragraph are of type (36). In order to determine their right-hand sides, we have to find the velocity components u and w, as well as space-derivatives of x and z.

3.2. The space-differentiating procedure

The derivatives of all quantities with respect to the generalized coordinates are calculated here using quadratic splines. This implies that some function f (C) defined on [C0, CM] is first approximated by a quadratic spline s(C):

f (C) « s(C) = a + bC + cC2, (42)

then the derivative of f is evaluated analytically:

f'(C) « s'(C) = b + 2 cC.

The spline coefficients a, b and c are piecewise constant: let ai-i/2, bi-i/2 and ci-i/2 (i = 1, 2,..., M) be their values within intervals (Ci-1, Ci) of some appropriate subdivision of [C0, CM]. To determine 3 M spline coefficients, we state the following requirements:

(i) s must be continuous:

ai-1/2 + bi-1/2 Ci + Ci—1/2 Ci2 = ai+1/2 + bi+1/2 Ci + Ci+1/2 Ci2 (M - 1 conditions );

(ii) s' must be continuous:

bi—1/2 + 2 Ci—1/2 Ci = bi+1/2 + 2 Ci+1/2 Ci (M - 1 conditions);

(iii) s must be equal to f in the middle points of intervals (Ci—1, Ci):

ai—1/2 + bi—1/2 Ci—1/2 + Ci—1/2 Ci2—1/2 = f (Ci—1/2) ( M conditions );

(iv) s must be equal to f in the end points of the whole interval [C0, CM]:

a1/2 + b1/2Co + C1/2Co = f(Co), aM—1/2 + bM—1/2Cm + cm—1/2Cm = f(Cm) ( 2 conditions).

The set of collocation points (in which the spline coincides with the function) is shown in figure 1 (a). Instead of f, its first or second derivatives may be posed in the end points, with the two last conditions being properly rewritten (but it is not necessary here). The resultant linear-equation system for one of three spline coefficients has a tridiagonal matrix which can be easily converted (by the Thomas algorithm). Then other two coefficients are calculated explicitly.

The use of quadratic splines seems to be most reasonable because they provide high enough (second) order of approximation and in the end points need the values of function only (in contrast to cubic splines for which both the function and its first or second derivative are necessary in the ends of interpolation interval).

3.3. The solution of the elliptic problem

To find the velocity components u and w, we have the system of interior equations (2), (3) with boundary conditions (13), (15), (17), (19) and corner conditions (24), (25), (27), (28), (30), (31), (33), (34). The coefficients J, a, ft, 7 involved in the system are calculated explicitly, according to (4), (5) and using the space-differentiating procedure described above. The equations (19), (31) and (34) should be first discretized in time, according to the time-stepping procedure.

Since the equation (2) is elliptic with respect to potential we can employ the iterative method of fractional steps [13]. More specifically, the stabilizing-corrector scheme [14, 13] adapted to equations of divergent type [15] is applied:

4>k _ k+1/2

= U

J

n

+wk,

^k+1 _ ^fc+l/2

(

= wf1

w

C '

Here u is an iteration parameter, the superscripts indicate iteration stages. To complete the system, the following equations are used:

u = + ^w)/a, w = (J^z + /^u)/y,

which are equivalent to (3) but seem more suitable for computation. Now the iterative procedure takes the form:

(i) the first (horizontal) half-step:

(u

k+l/2

+jwk,

uk+1/2 - (J/a) ri/2 = (£/a) wk;

(ii) the second (vertical) half-step:

k+1 - (WJ+1

(w

w^1 - (J/7) = (£/7) uk+1/2.

(43)

(44)

(45)

(46)

Both the systems (43) - (44) and (45) - (46) are one-dimensional: they are differential with respect to one generalized coordinate whereas another coordinate is involved only as a parameter (the derivatives in the right-hand sides of equations are known from the previous iteration stage). This parameter being discretized, we have the set of one-dimensional problems at every half-step. All these problems are calculated with use of quadratic splines. Both the unknown functions, namely ^ and u (let us consider the horizontal half-step for example), are approximated by quadratic splines of kind (42). Two equations are satisfied in each node that yields 2 (M + 2) conditions for determining 6 M spline coefficients. The lacking 4 (M — 1) conditions are provided by the continuity requirements for the splines and their first derivatives. The operation order at the vertical half-step is quite similar.

The two-dimensional version of the set of collocation points is shown in figure 1 (b).

k

Figure 1. . The one-dimensional (a) and two-dimensional (b) sets of collocation points.

Table 1.

Nodes Horizontal half-step Vertical half-step

1 (43) + (44) (45) + (46)

2 (43) + (15) or (17) (45) + (46)

3 (43) + (44) (45) + (13)

4 (43) + (44) (46) + (19)d

5 (43) + (24) or (27) (45) + (25) or (28)

6 (43) + (30) or (33) (46) + (31)d or (34)d

All the collocation points are subdivided into six types, depending on equations considered to be valid in them. These equations are listed in table 1, where the superscript d indicates that the evolutionary equations (19), (31) and (34) should be first discretized in time. The iteration interrupts when the following conditions become achieved:

II k+1/2 . k || ^ II k+1/2 . k+1 || „

|| u + w]k || < e, || u + w>k+1 II < e,

with some sufficiently small e > 0. The norm means a maximal absolute value, with the maximum being taken over the nodes 1-3, 5 (in which the equations (43) and (45) are both valid).

4. Results and discussion

To test the numerical method developed, the following three problems have been calculated: on a liquid ellipse, on an overturning bore (TANH-wave) and on a disturbed solitary wave of large amplitude. Possessing their own advantages and shortcomings, these problems complete each other to make it possible to test the numerical method rather comprehensively.

4.1. A liquid ellipse

The solution of this problem describes the water motion in the absence of gravitation when all the fluid boundary is free and coincides with an ellipse:

S(x,z,t) = Q2 + (az)2 - 1 = 0, (47)

where a = a(t) is a horizontal ellipse half-axis. This problem represents a saldom example of free-surface problems possessing an exact solution in nonlinear formulation. More precisely, the system of partial differential equations is reduced to a single ordinary differential equation governing the evolution of a(t) [16] (dots above stand for time derivatives):

a2

a = c . , a(0) = ao

VTTa4 ^ 0

with c being an arbitrary constant.

It is easy to check that the following potential (the tilde indicates an exact solution)

(48)

<§(x,z,i)

1 a

2 a

(x2 - z2)

(49)

is harmonic inside the fluid domain and satisfies the kinematic condition

St + Sx + Sz = 0,

as well as the dynamic condition

*t + 1($x + *2) = c (t),

on the free boundary (47). The latter condition is Bernoulli's equation with zero gravitation (g = 0) but non-zero right-hand function, namely:

C(t) =

•2 / 2 \ 2 a2 / ca2

1 + a4 V 1 + a

The fluid-particle trajectories can be determined by integrating (7) with $(x,z,t) given by (49)

X(f, Z, t) = Xo(f, Z) —, z(f, Z, t) = zo(f, Z) . (50)

-0 -(t)

The problem is adapted to a quadrangular fluid domain by setting:

L(x,z,t) = x, R(x,z,t) = x — X (t), B(x,z,t) = z.

This means that the left wall and bottom coincide with the ellipse symmetry axes. The right wall is vertical and movable (a piston): the law of its motion is defined according to the first of (50):

X (t) = Xo ^,

-o

with some constant X0 (0 < X0 < -0).

The initial value of potential follows directly from (49), (48):

<M0 = $o(XoN (f), Z0N (f)),

1 ai / 2 2\ • a0_

$0(x, z) = - — (x — z ), ai = a(0) = c

2 ao ' v/TT

o

x0N(£) = xo(£, (N), zoN(£) = z0(£ (N) = — \ 1 - xoN(£)

- o - o

The initial variable transformation is chosen to be rather simple:

Xo(e-Z) = Xo M zoie.c) = zoN<e) N = 1 — (X s)2i

The following values of flow parameters

-0 = 1, -1 = — 1 ( ^ c = — v^), Xo = 0.5

are used in calculations carried out on the grid M x N with (i) M = N = 10, (ii) M = N = 20, (iii) M = N = 30. The time step t is chosen to amount to 0.1/M.

Table 2.

t II x || || x — x || / || x || x 104 " * " || z — z || / || z || x 104

10 x 10 20 x 20 30 x 30 10 x 10 20 x 20 30 x 30

0 0.500 0.00 0.00 0.00 1.000 0.00 0.00 0.00

1 0.220 0.64 0.18 0.08 2.276 3.24 0.82 0.37

2 0.136 1.57 0.53 0.26 3.679 5.14 1.30 0.58

3 0.098 2.62 0.96 0.50 5.091 6.44 1.62 0.74

4 0.077 3.87 1.53 0.79 6.505 7.44 1.88 0.86

5 0.063 5.42 2.25 1.08 7.919 8.25 2.09 0.96

6 0.054 7.34 3.09 1.32 9.333 8.96 2.27 1.04

7 0.046 9.68 4.02 1.77 10.747 9.59 2.42 1.12

The results of calculations are presented in table 2. Here the norm of any function means its maximal absolute value (in space) so that

II ~ II II II a Xo || ~ | | | 1

ao ao a a

The problem has been calculated up to the moment t = 7, when the vertical size of fluid domain becomes more than 10 times greater (and, correspondingly, the horizontal size more than 10 times smaller) so that the vertical-to-horizontal size ratio becomes more than 100 times greater compared to its initial value. One can see that the relative errors increase with time but do not exceed 0.1% in the worst case. The relative errors decrease with increasing the number M x N of grid sells. The actual order of decreasing is lower than the theoretical (second) one, but closer to the second than to the first.

4.2. An overturning bore (TANH-wave)

To demonstrate the numerical method capability of treating overturning waves, we have chosen the so-called TANH-wave [17]. This wave can be generated from the following initial free-surface configuration

Z0N (f) = n(x0N (f))

1 1 / r — b

n(x) = |U| + - U2, U(r) = - c 1 + tanh r b

4 ' w 2 V A

and surface value of potential:

<Mf) = $o(roN (0),

('x 1 ( r - b $o(r) = $o(0) + U(y) dy =- c r + A ln cosh—-— o 2 A

Here roN (f) is some appropriate initial distribution of grid nodes on the free surface.

Far enough from the centre of TANH-wave (r = b) the water motion is close to a homogeneous horizontal flow with velocity equal to 0 at the left and to c at the right. Hence, the following vertical walls can be placed:

L(r, z, t) = r = 0, R(r, z, t) = r — (l + ct) = 0.

From below the water is bounded by a horizontal bottom:

B (r,z,t) = z + h = 0.

Figure 2. . t = 3.5; 4.1;

An 4.7;

overturning bore 5.3 (from right to

(TANH-wave): a — the free-surface configurations at moments left); b — the grid node locations at moment t = 5.0.

The following values of flow parameters

g = 1, h = 1, l = 30, c = - 1.7, b =15, A = 2, and grid parameters

M =150, N = 10, t = 01

M

are used in calculations.

The results of calculations are presented in figure 2. The free-surface configurations at moments t = 3.5; 4.1; 4.7; 5.3 (a) and the grid node locations at moment t = 5.0 (b) are exhibited. Squares indicate the actual positions of fluid particles (grid nodes). By means of dotted curves, the fluid-particle trajectories (a) and the coordinate lines C = const (b) are also shown for the sake of better visualization.

One can see that by the final moment t = 5.3 (just after which the computations fail) the TANH-wave developes to a pronounced "plunging breaker" (see [18] for the breaking-wave classification). The wave profile steepens and, after the surface tangent becomes vertical at some point, overturns. A prominent jet of water forms and then falls down in front of the wave. Such fall usually causes a large splash, but here it is absent because the jet "does not know" about its meeting with still water. For all "life period" of the wave, its profile remains smooth (with no artificial smoothing). The computational-grid cells stretch essentially in vertical direction. Such strong grid deformation seems to be the reason why the calculations become impossible to be continued.

4.3. A disturbed solitary wave of large amplitude

A solitary wave is known to be some specific solution of the potential model. This solution is essentially nonlinear and dispersive. In addition, it is determined by a single parameter (by a wave amplitude usually). All this makes a solitary wave very attractive for testing numerical algorithms.

A free surface can be represented in the form:

z = n(x, t) = a f (X), X(x, t) = X - (bA+ ct), (51)

with a, b, c and A being the amplitude, initial phase, speed and length of a solitary wave. The parameter b is unessential, the remaining three are related with each other. These relations, as well as the kind of function f (X), are specified below, now we only mention the necessary properties of this function: (i) f (0) = 1,

(ii) f (-X ) = f (X),

(iii) Xf'(X) < 0,

(iv) f (X) ^ 0 as X ^ (with some ^ > 0) .

The property (iv) indicates that water is disturbed only in the vicinity of a solitary-wave crest but actually undisturbed far enough from it. So the fluid domain can be bounded by fixed vertical walls:

L(x, z, t) = x = 0, R(x, z, t) = x — l = 0.

These walls do not affect a solitary wave during some period which is proportional to b or l — b, depending on the direction of wave propagation. From below the water is bounded by a level bottom:

B (x,z,t) = z + h = 0. The initial free-surface configuration follows directly from (51):

/ i/-\ n\ !■ I X0N ) — b

z0N(C) = n(x0N^ 0) = af

A

where x0N (£) is some appropriate initial distribution of grid nodes on the free surface. As for the initial surface value of potential, it can be obtained in the following manner. Based on the obvious consequence from (51), namely

Xt + cXx = 0,

one can derive the following (steady) free-surface condition:

w = — cz? (Z = Zn ). (52)

This condition makes it possible to find the initial value of potential (in the whole domain) as a solution of the elliptic problem (2), (3), (13), (15), (17), (52), (24), (25), (27), (28). It should be emphasized that the steady condition (52) is used only at instant t = 0, then it replaced by the usual unsteady condition (19).

Within the potential model, the function f (X) can not be represented in closed form. Available are only some approximate expressions for it, the simplest of which is (see [19])

f (X ) = sech2 X. (53)

The corresponding approximate expression for the wave length A looks as

A\2 4 h

V 3 a- (54)

The speed c can be obtained from the integral (exact) properties of solitary waves (see [20]):

— = 1 + 3 ah (55)

gh = 1 + 2 hh , (55)

/OO rtt

f (X) dX, /2 = f 2(x) dX.

-O J —O

Since the relations (53), (54) and, as a consequence, (55) are not exact within the potential model, the initial wave defined by these relations should be considered as "a disturbed solitary wave".

The following values of flow parameters

g = 1, h = 1, l = 50, b =35,

and grid parameters

M = 250, N = 10, t =

M

34 X

Figure 3. . Solitary waves with amplitudes a = 0.6; 0.7; 0.8; 0.9; 1.0 (from top to bottom).

are used in calculations. The amplitude - has been varied in the neighbourhood of value - = 0.833h which is maximal for solitary waves (see [21]). The speed c and length A have been evaluated according to the relations (55) and (54) simplified by taking into account that g = h =1 and I1 = 2, 12 = 4/3 for the function f (X) given by (53):

/- > 2

c = — v 1 + -, A

\f3a

The sign " - " in the former relation indicates that the solitary wave is prescribed to travell from right to left.

The results of calculations are presented in figure 3. Shown are the profiles of solitary waves with amplitudes a = 0.6; 0.7; 0.8; 0.9; 1.0 (from top to bottom). The figure has qualitative rather than quantitative character, let us therefore not concretize the time increments with which the profiles are plotted. Notice only that the latest time moments (the left profiles) are closest to the moments of computation failure.

One can see that only the lowest from the considered waves remains solitary, all the others break. The larger is the wave amplitude, the more distinctive is the overturning phase. That is, the amplitude being increased, the wave transforms gradually from "the spilling breaker" (a = 0.7) to "the plunging breaker" (a > 0.9), according to the commonly accepted classification (see [18]).

The disturbances of solitary waves appear crucial, as the present numerical results demonstrate, for the waves with amplitudes a > 0.7. This sligthly differs from the results of Tanaka [22], according to which a solitary wave becomes unstable at value 0.78 of its amplitude. Probably, the stability radius of a solitary wave with a = 0.7 is too small so that the disturbance applied has brought the wave out this radius.

The results of test calculations allow us to conclude that the numerical method developed is capable of treating, with high enough accuracy, strongly nonlinear waves, including overturning ones (from plunging to spilling breakers).

This work was partially supported by the Sibirian Branch of Russian Academy of Sciences (SD RAS) within the Integrate Project No. 43, and by the Royal Society of London, together with the SD RAS, within the Joint English-Russian Project "Non-linear Plane Problem on Water Impact".

References

[1] LONGUET-HlGGINS M.S., COKELET E. D. The deformation of steep surface waves on water. I. A numerical method of computation. Proc. R. Soc. Lond. A, 350, 1976, 1-26.

[2] VlNJE T., BREVIG P. Numerical simulation of breaking waves. Adv. Water Resources, 4, 1981, 77-82.

[3] Baker G. R., Meiron D. I., Orszag S. A. Generalized vortex methods for free-surface flow problems. J. Fluid Mech., 123, 1982, 477-501.

[4] Dold J.W., PEREGRINE D. H. An efficient boundary-integral method for steep unsteady water waves. In "Numerical Methods for Fluid Dynamics II" (ed. K. W. Morton, M.J. Baines). Clarendon Press, Oxford, 1986, 671-679.

[5] COINTE R. Nonlinear simulation of transient free surface flows. In "Proc. 5th Intl Conf. Numer. Ship Hydro." Hiroshima, 1989, 239-250.

[6] Grilli S.T., Skourup J., Svendsen I. A. An efficient boundary element method for nonlinear water waves. Eng. Anal. Bound. Elem., 6, 1989, 97-107.

[7] Dold J. W. An efficient surface-integral algorithm applied to unsteady gravity waves. J. Comp. Phys, 103, 1992, 90-115.

[8] HAUSSLING H.J. Solution of nonlinear water wave problems using boundary-fitted coordinate systems. In "Numerical Grid Generation" (ed. J. F. Thompson). North-Holland, N. Y., 1982, 385-407.

[9] Telste J. G. Calculation of fluid motion resulting from large-amplitude forced heave motion of a two-dimensional cylinder in a free surface. In "Proc. 4th Intl Conf. Numer. Ship Hydro." Washington, 1985, 81-93.

[10] Yeung R. W., Vaidhyanathan M. Non-linear interaction of water waves with submerged obstacles. Int. J. Numer. Meth. Fluids, 14, 1992, 1111-1130.

[11] Wu G.X., EATOCK Taylor R. Finite element analysis of two-dimensional non-linear transient water waves. Appl. Ocean Res., 16, 1994, 363-372.

[12] Wu G. X., EATOCK Taylor R. Time stepping solutions of the two-dimensional nonlinear wave radiation problem. Ocean Eng., 22, 1995, 785-798.

[13] YANENKO N. N. The Method of Fractional Steps for Solving Multi-Dimensional Problems of Mathematical Physics. Nauka, Novosibirsk, 1967 (in Russian; transl. to English: Springer-Verlag, Berlin, 1971).

[14] DOUGLAS J., Rachford H. H. On the numerical solution of heat conduction problems in two and three space variables. Trans. Amer. Math. Soc., 82, 1956, 421-439.

[15] PROTOPOPOV B. E. Numerical simulation of wave motions in ideal fluid of a finite depth. Trans. Soc. Naval Archit. Korea, 32, 1995, 58-69.

[16] OVSYANNIKOV L. V. General equations and examples. In "The Problem on Unsteady Motion of Fluid with a Free Boundary". Nauka, Novosibirsk, 1967, 5-75 (in Russian).

[17] COOKER M. J. The Interaction of Steep Water Waves and Coastal Structures (Ph. D. Thesis). University of Bristol, 1990.

[18] Peregrine D. H. Breaking waves on beaches. Ann. Rev. Fluid Mech., 15, 1983, 149-178.

[19] Miles J.W. Solitary waves. Ibid, 12, 1980, 11-43.

[20] LONGUET-HlGGINS M. S. On the mass, momentum, energy and circulation of a solitary wave. Proc. R. Soc. Lond. A, 337, 1974, 1-13.

[21] EVANS W. A. B., Ford M. J. An exact integral equation for solitary waves (with new numerical results for some 'internal' properties). Ibid, 452, 1996, 373-390.

[22] TANAKA M. The stability of solitary waves. Phys. Fluids, 29, 1986, 650-655.

Received for publication October 9, 1997