Numerical simulation of non-linear regular and focused waves in an infinite water-depth

ARTICLE IN PRESS

Ocean Engineering 35 (2008) 887–899 www.elsevier.com/locate/oceaneng

Numerical simulation of non-linear regular and focused waves in an infinite water-depth D.Z. Ninga,b, B. Tenga, R. Eatock Taylorb,, J. Zangc a

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116023, China b Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, UK c Department of Architecture and Civil Engineering, University of Bath, Bath BA2 7AY, UK Received 4 January 2007; accepted 16 January 2008 Available online 31 January 2008

Abstract Inviscid three-dimensional free surface wave motions are simulated using a novel quadratic higher order boundary element model (HOBEM) based on potential theory for irrotational, incompressible fluid flow in an infinite water-depth. The free surface boundary conditions are fully non-linear. Based on the use of images, a channel Green function is developed and applied to the present model so that two lateral surfaces of an infinite-depth wave tank can be excluded from the calculation domain. In order to generate incident waves and dissipate outgoing waves, a non-reflective wave generator, composed of a series of vertically aligned point sources in the computational domain, is used in conjunction with upstream and downstream damping layers. Numerical experiments are carried out, with linear and fully non-linear, regular and focused waves. It can be seen from the results that the present approach is effective in generating a specified wave profile in an infinite water-depth without reflection at the open boundaries, and fully non-linear numerical simulations compare well with theoretical solutions. The present numerical technique is aimed at efficient modelling of the non-linear wave interactions with ocean structures in deep water. r 2008 Elsevier Ltd. All rights reserved. Keywords: Numerical wave tank; Infinite water-depth; HOBEM; Focused wave groups; Fully non-linear

1. Introduction The use of the numerical wave tank has been developed during the last two decades as an effective and efficient tool for modelling various water wave problems. For example, Romate (1989), Xu (1992), Kim and Kim (1997) and Ning and Teng (2007) employed higher order boundary element methods to simulate fully non-linear waves. Brandini and Grilli (2001) and Fochesato et al. (2007) considered the numerical simulation of rogue waves using a 16-node higher order boundary element model (HOBEM). Ferrant (1998) and Ferrant et al. (2003) used similar methods for the fully non-linear diffraction problem in regular and irregular waves. Celebi (2001), Scorpio et al. (1996) and Beck (1994) adopted the desingularized boundary integral equation method to calculate fully non-linear wave loads Corresponding author. Tel.: +44 1865 273144; fax: +44 1865 273010.

E-mail address: [email protected] (R. Eatock Taylor). 0029-8018/$ - see front matter r 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.oceaneng.2008.01.015

on objects; Wu et al. (1998), Hu et al. (2002) and Wang and Khoo (2005) used the finite element method to study fully non-linear water wave problems. Xu et al. (1993) used the Green–Naghdi theory for the numerical simulation of irregular non-linear waves, and Turnbull et al. (2003) adopted a s-transformed finite element method to study 2D fully non-linear water wave problems. A more detailed review of some of this work has been given by Kim et al. (1999). All of these studies are concerned with the case of finite water-depth. In certain conditions, it could be computationally more efficient, and physically realistic (such as Ferrant, 1990; Lee et al., 1994), to model problems as being in infinite water-depth. This is the motivation for the present work. In this paper, an efficient HOBEM is developed to simulate a non-linear numerical wave tank in an infinite water-depth without a structure. The formulation is 3D as this study is a precursor to an investigation of the corresponding fully non-linear wave–body interaction problem in deep water.

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

888

In order to develop a numerical wave tank using boundary elements in an infinite water-depth domain, two key problems must be dealt with. One is the realization of an appropriate means of wave generation; and the other is the exclusion of the infinite-depth lateral surfaces of the wave tank. In this paper, an efficient 3D HOBEM is developed to resolve these issues. The approach is an extension of the work of Ning and Teng (2007), which dealt with non-linear regular and irregular waves in finite waterdepth. A channel Green function in infinite water-depth is used here based on image theory, so that the two lateral surfaces are excluded from the computational domain and the proposed problem is transformed from an infinite domain to a finite domain. A wavemaker based on submerged sources is used here to generate arbitrary input waves, with a starting ramp applied in order to reduce the effect of transient disturbances. At the ends of the tank, artificial absorbing layers are employed on the free surface to act as beaches. The boundary integral equation is formulated with quadratic elements and the resulting matrix equation system is then solved at each time step. The instantaneous free surface is updated at each time step by using the fourth-order Runge–Kutta time integration scheme with the semi-mixed Eulerian–Lagrangian time stepping procedure. Calculations are made concerning the generation and propagation of fully non-linear regular waves and focused wave groups in deep water. The numerical results show that the proposed scheme possesses high accuracy and good numerical stability. The following sections first introduce the governing equation and boundary conditions, and the higher order boundary integral equation based on the image Green function, and then give results from a series of numerical experiments.

2. Governing equation and boundary conditions A Cartesian coordinate system is defined with the origin in the plane of the undisturbed free surface z ¼ 0, with the z-axis positive upwards, as shown in Fig. 1. The corresponding boundary value problem defining the fluid motion is set up as follows. Let t denote time and Z the free surface elevation above the still water level. The waterdepth is infinite, and the flow is 3D. The free surface is Sf and the solid boundaries of the fluid domain (the lateral surfaces) are denoted SB. The waves are generated by z Damping zone

Damping zone

xb1

xb2

o

Lb h

sources

h∞

Lb

~ v ¼ rf,

(1)

where ~ v ¼ ðu; v; wÞ, r ¼ ðq=qx; q=qy; q=qzÞ, and u, v, w are the velocities in the x-, y- and z-directions, respectively. Generally, the velocity potential satisfies the Laplace equation. However, due to the presence of the source distribution, the velocity potential satisfies the Poisson equation in this case (Brorsen and Larsen, 1987): r 2 f ¼ q , 2

(2) 2

2

2

2

2

2

where r ¼ q =qx þ q =qy þ q =qz is the 3D Laplacian operator, and q ¼ q ðx; y; zÞ is the pulsating volume flux density of the internal source distribution described below. In other words, the volume flux through an infinitesimal area surrounding a volume dx dy dz of the distributed source distribution is q dx dy dz. This approach for generating incident waves differs from that usually adopted, as, for example, by Ning and Teng (2007). It has the advantages of generating arbitrary waves and allowing dissipation of the waves reflected by objects back to the input boundary. On the instantaneous free surface, the fully non-linear kinematic and dynamic boundary conditions can be derived from the material derivatives and Bernoulli’s equation. In the present work, the semi-mixed Eulerian– Lagrangian approach is used, in which free surface nodes are allowed to move in the vertical direction with nodal velocity (~ v). The resulting equations on the free surface Sf can be written as follows in the semi-Lagrangian frame: dZ qf qZ qf qZ qf ¼  þ ; dt qx qx qy qy qz

on Sf

(3)

and df 1 dZ df ¼  rf  rf  gZ þ ; on S f , (4) dt 2 dt dz where g is the gravitational acceleration and the time derivative is d q ¼ þ~ vr dt qt with the velocity vector defined by ~ v ¼ ð0; 0; dZ=dtÞ. At the lateral walls of the tank, the zero normal velocity condition is imposed: qf ¼ 0; on S B , (5) qn where n is the outward normal direction of the solid boundary. Since the proposed problem is solved in the time domain, an initial condition must also be imposed. Here, we specify vn ¼

h∞

Fig. 1. Definition sketch.

x

controlling the volume flux density of a vertical source distribution inside the model boundary. It is assumed that the fluid is incompressible, inviscid, and the flow irrotational. The fluid motion can therefore be described by a velocity potential f, related to the velocity by

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

that the water surface is initially calm, so that j ¼ Z ¼ 0;

tp0.

(6)

At each subsequent time step, a Dirichlet condition is applied on the free surface, using the dynamic boundary condition to update the value of the potential from the previous time step. 3. Higher order boundary element method 3.1. Boundary integral equation By applying Green’s second identity to the fluid domain O, the boundary value problem presented in the previous section can be converted in the usual manner into the following boundary integral equation (e.g., Brebbia and Walker, 1980; Ligget and Liu, 1983)  Z  qGðp; qÞ qfðqÞ aðpÞfðpÞ ¼  Gðp; qÞ fðqÞ dS qn qn S Z þ q Gðp; qÞ dO, (7) O

where p ¼ (x0, y0, z0) and q ¼ (x, y, z) are source and field points, respectively, and a(p) is the solid angle. G is the Green function. It would be difficult to solve the present problem if the simple 1/r Green function were used because of the infinite computational domain. Here, a channel Green function is adopted, which can be obtained by the superposition of infinite images reflected about the two lateral walls. They can therefore be excluded from the computational domain. The Green function is then written as (Newman, 1992): " 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Gðp; qÞ ¼  2 4p X þ Y 2 þ Z2 0 1 X 1 B þ @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n¼1 X 2 þ ðY þ 2nBÞ2 þ Z 2 13 1 1 C7 þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  A5, nB X 2 þ ðY  2nBÞ2 þ Z 2

(8)

where X ¼ xx0, Y ¼ yy0, Z ¼ zz0 and B is the width of the tank. To ensure the convergence of the Green function, a factor 1/nB is subtracted from each term (Newman, 1992). The terms in the summation in Eq. (8) represent reflections of the sources in the lateral walls. A similar technique was used to solve linear wave diffraction/radiation problems in the frequency domain (Teng and Eatock Taylor, 1995). Following Brorsen and Larsen (1987), the incident wave is specified by writing the flux density in Eqs. (2) and (7) as follows: q ¼ 2vdðx  xs Þ.

(9)

Here, v=v(x, y, z, t) is the horizontal fluid speed corresponding to the wave to be generated, xs is the horizontal

889

position of the vertical source and d(xxs) is the Dirac delta function. 3.2. Discretization The boundary surface is discretized by higher order elements based on quadratic shape functions, as discussed by Ning and Teng (2007). Within the boundary elements, physical variables are also interpolated by the same shape functions, i.e., the elements are isoparametric. The resulting approximations are ½x; y; z ¼

K X

hk ðx; BÞ½xk ; yk ; zk ;

fðx; BÞ ¼

k¼1

K X

hk ðx; BÞfk ,

k¼1

  qf qf k ¼ h ðx; BÞ , qn qn k k¼1 K X

(10)

where (x, B) represent the local intrinsic coordinates; [xk, yk, zk], fk, ðqf=qnÞk and hk are the coordinates, potential, normal derivative of the potential and the shape function corresponding to the kth node; and K is the number of nodes in the element. The shape functions can be found in standard finite element texts such as Zienkiewicz and Morgan (1983). By substituting Eqs. (9) and (10) into Eq. (7), the integral equation is obtained as follows: aðpÞfðpÞ 

N e2 Z X i¼1

i¼1

þ

Z

1

XZ N e2

¼ 

1

1

K X

hk ðx; BÞfk

1 k¼1

Z

1

Gðp; qðx; BÞÞ

1

XZ

1

N e1

1

i¼1

1

1

Z

qGðp; qðx; BÞÞ jJðx; BÞj dx dB qn

qfðqðx; BÞÞ jJðx; BÞj dx dB qn

1

2vGðp; qðx; BÞjJðx; BÞj dx dB,

(11)

1

where J(x, B) represents the Jacobian matrix relating the global coordinates and the local intrinsic coordinates in the ith element. N e1 and N e2 are the numbers of discretized elements on the source surface and the free surface, respectively. The above single- and double-layer integrals are calculated numerically using the 16-point Gauss rule if p6¼q. If p ¼ q, a polar coordinate transformation (Li et al., 1985) is employed to evaluate the single-layer singular integrals, and a direct method of subdividing the singular integral into two parts, one including the singularity and the other being regular, is applied to estimate the doublelayer singular integrals. The term including the singularity is dealt with by accumulating all the neighbouring singular integrals surrounding the integration point, as explained by Teng, et al. (2006). It can be seen from the second term of the right hand side in Eq. (11) that wave sources are distributed on all the Gauss points of each source surface element. The final system of equations is obtained in a standard way by assuming that the discretized equations are satisfied exactly at a defined set of collocation points. This leads to

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

890

the matrix equation: A x ¼ b,

(12)

where x is a vector of unknown normal velocities on the free surface; A is the influence coefficient matrix and b is the vector obtained from the right hand side of Eq. (11). These equations are reformulated at each time step by a fourth-order Runge–Kutta integration scheme, as the nodes are moved vertically and the grid is updated by stretching. No regridding as such is required. The methodology for obtaining the velocity components of the nodes, the updated position of the free surface, and the updated velocity potential is the same as described by Ning and Teng (2007). The difference is that the vertical positions of the sources are also updated following the changes of the free surface position, the lowest point of the source surface being fixed. The vertical extent of the source surface is discussed below. At the start of the computation, a cosine ramp function is applied over a time 2T, where T is the period of the incident regular wave, or the minimum period of wave components in a focused wave group. At each end of the domain (behind the wave-making sources and far downstream), a numerical beach is created by means of the artificial damping layer. Additional dissipative terms are included in the non-linear dynamic and kinematic free surface boundary conditions, again as described by Ning and Teng (2007). 3.3. Smoothing When steep waves are simulated it is found that smoothing on the free surface is required in order to prevent so-called saw-tooth instabilities. This involves assigning new values to the nodal values of Z and j based on a weighted average of the initially calculated nodal values. The following five-node smoothing function is used: Zi ¼

1 ðZi2 þ 4Zi1 þ 10Zi þ 4Ziþ1  Ziþ2 Þ 16

(13)

and 1 ðfi2 þ 4fi1 þ 10fi þ 4fiþ1  fiþ2 Þ. (14) 16 This scheme was adopted by Longuet-Higgins and Cokelet (1976), and similar approaches have been used by several other investigators to overcome the instabilities that can arise in wave problems. In our calculations, smoothing is used at different times according to the different wave slopes arising.

fi ¼

4. Numerical results and discussions We now illustrate the method described by undertaking several numerical simulations of linear and non-linear, regular and focused waves in the proposed numerical wave tank, using an analytical input wave velocity at the source surface. First, some numerical experiments were performed

to check the convergence of the Green function. These showed that converged solutions can be obtained from Eq. (11) when the infinite sums in Eq. (8) are truncated at 20 terms, for the width of tank, B, used in the present paper. 4.1. Simulation of linear and non-linear regular waves Progressive free surface waves are generated in the present numerical tank by imposing the horizontal water particle velocity components on the sources distributed vertically in the middle of the tank (see Fig. 1). For linear waves, the horizontal velocity component required in Eq. (11) is calculated from the usual linear equation: vðx; y; z; tÞ ¼ Ao ekz cosðkx  otÞ.

(15)

The corresponding free surface elevation is ZI ðx; y; z; tÞ ¼ A cosðkx  otÞ,

(16)

where k is the wave number, A the wave amplitude, o the wave frequency, and k and o satisfy the deep water dispersion equation: o2 ¼ gk.

(17)

For non-linear waves, the incident horizontal velocity is here determined from the fifth-order Stokes analytical velocity potential for deep water (Fenton, 1985):   1 1 37 vðx; y; z; tÞ ¼ k  k3  k5 ekz cosðkx  otÞ k 2 24 1 þ k4 e2kz cos 2ðkx  otÞ þ k5 e3kz 4 rffiffiffi g  cos 3ðkx  otÞ . (18) k The corresponding free surface elevation is given by     1 3 422 5 1 2 1 4   3   cosðkx  otÞ þ  þ  ZI ¼ k 8 384 2 3   3 3 99 5  þ  cos 3ðkx  otÞ  cos 2ðkx  otÞ þ 8 128  1 4 125 5  cos 5ðkx  otÞ , (19) þ  cos 4ðkx  otÞ þ 3 384 where the wave slope e ¼ kA satisfies: 1þ

2 4 o þ ¼ pffiffiffiffiffiffi . 2 8 gk

(20)

From Eqs. (15) and (18), it can be seen that both linear and fully non-linear velocities decrease exponentially with increasing water-depth. That is to say, sources contribute little to the wave generation when the depth of the source exceeds some value, and the corresponding components then tend to zero. That value is here taken as 1.5 times the wavelength, following some numerical experiments. The first example concerns the propagation of a regular wave with frequency o=3.14 rad/s corresponding to a 2 s wave in a flume. As indicated above, we are concerned with

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

1.5 coarse mesh

intermediate mesh

fine mesh

1

η/A

0.5 0 -0.5 -1 -1.5 6

7

8

9

10

11 t(s)

12

13

14

15

16

Fig. 2. Time history of wave elevation at (l/2, B/2) for different element sizes.

1.5 linear theory

Present

1

η/A

0.5 0 -0.5 -1 -1.5 0

4

8

12

t/T Fig. 3. Time history of wave elevation at (1.5l, B/2).

1.5 10T

11T

linear theory

1 0.5 η/A

the problem in deep water. The corresponding linear and fully non-linear (fifth-order) wave number k and wavelength l can be obtained from Eqs. (17) and (20), respectively. The computational domain is defined as L (length)=6l(1.5lpxp4.5l) and B (width)=0.2l(0.0pyp0.2l), respectively. To absorb the propagating wave energy, damping layers of length 1.5l are placed at both ends of the domain, these parameters being selected to ensure that all outgoing waves can be dissipated, as demonstrated by Teng and Ning (2006). Numerical experiments were initially carried out to assess the convergence of the adaptive spatial and temporal discretization for a wave slope kA ¼ 0.09. Fig. 2 depicts the free surface time history at the position (l/2, B/2) obtained using a fine mesh (24 quadratic elements per wavelength, i.e., Dx ¼ l/24); an intermediate mesh (16 elements per wavelength, i.e., Dx ¼ l/16); and a coarse mesh (eight elements per wavelength, i.e., Dx ¼ l/8). In all cases here, two elements were used across the width of the numerical wave tank and 30 elements are distributed on the source surface. The results obtained using the fine and intermediate meshes are almost identical, indicating that mesh convergence was achieved using the intermediate mesh (Dx ¼ l/16). Similar tests for different time steps (T/20, T/40 and T/60) using the intermediate mesh indicate that Dt ¼ T/40 is sufficient. Fig. 3 compares the wave elevation at location (1.5l, B/2.0), plotted against dimensionless time in the case of linear slope kA ¼ 0.02, with the linear theoretical solution, the wave elevation being normalized by the amplitude of the incident wave. It can be seen from the figure that the initial wave amplitude growth is smooth due to the ramp function, and the computed wave agrees very well with the linear theory when the steady state has been reached. Surface elevations along the centre line of the wave tank at 10T, and 11T are plotted in Fig. 4. The effects of the damping layers are indicated by the decreasing wave height over the damping zones at each end. The wave fields at these two times are almost identical, suggesting that the steady state has been reached at t/T ¼ 10. The wave crests and troughs at these times are very close to the predictions from linear theory. For this case, no smoothing was used.

891

0 -0.5 -1 -1.5 -10

-5

0

5

10 x(m)

15

20

25

30

Fig. 4. Wave elevations along the plane of symmetry.

By making use of Eq. (18) as the incident velocity on the input boundary and fully non-linear boundary conditions satisfied on the free surface, steeper wave non-linear propagation in infinite water-depth is simulated for several cases. The time histories of the simulated fully non-linear wave elevations, at locations (0.5l, B/2) for kA ¼ 0.21 and (1.0l, B/2) for kA ¼ 0.31, are presented in Fig. 5 and compared with the fifth-order theoretical input waves. In each case the fundamental wave frequency is o ¼ 3.14 rad/s, as above. The numerical results based on smoothing are in good agreement with the fifth-order theory. They show that the wave amplitudes are almost unchanged after the initial development of the waves, and the non-linear features are exhibited very clearly, such as higher crests and smaller troughs. The non-linear wave features become more evident with increasing wave slope kA. The results also confirm that the phase error is very small, as found for the linear waves. It is of interest to investigate the performance of the model in representing different order non-linear wave components, such as linear, second- and third-order, etc. Taking Fig. 5(b) as an example, the contributions from the lowest order wave components have been calculated by use of the FFT

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

1.5

1.5

1

1

0.5

0.5 η/A

η/A

892

0

0

-0.5

-0.5

-1

-1 -1.5

-1.5 0

2

4

6 t/T

8

10

12

Present

0

2

4

6 t/T

8

10

12

5th-order Stokes theory

Fig. 5. Time histories of surface elevation with different wave slope.

There are no parasitic transverse waves, which are known to be a potential source of difficulty, and this helps to confirm the successful implementation of the image Green function.

200 Present 5th-order Stokes

160

4.2. Simulation of linear and fully non-linear focused waves

80

Focused wave groups can be used to represent the profiles of extreme waves. In such a group, often referred to as a NewWave, all wave components in a spectrum focus simultaneously at a position in space. The representation has been studied experimentally and numerically by several investigators, such as Tromans et al. (1991), Taylor and Haagsma (1994), Baldock et al. (1996) and Borthwick et al. (2006). We now use the HOBEM discussed here, combined with source generation, to simulate linear and fully nonlinear focused wave groups in an infinite water-depth. For the linear focused wave group, the horizontal velocity component in Eq. (9) is obtained by superposition as

S(f)

120

40

0 0

0.5

1

1.5 f(Hz)

2

2.5

3

Fig. 6. Distribution of spectrum from different wave components.

(fast Fourier transformation), and they are plotted against the frequency f in Hz in Fig. 6. The different contributions can be seen clearly, up to the fourth-order wave. The percentages of the total arising from first- to fifth-order wave components are calculated as 81.26%, 13.64%, 3.53%, 1.13% and 0.44%, respectively. The surface elevations of the simulated wave along the centre line at t ¼ 10T and 13T for the case kA ¼ 0.21, and at t ¼ 10.5T and 12.5T for the case kA ¼ 0.31, are illustrated in Fig. 7(a) and (b), respectively. The corresponding fifth-order Stokes theoretical solutions are also plotted in these figures. We find that reflections from the upwave and downwave boundaries are not perceptible in these cases of steady wave propagation. Fig. 8 shows the comparison of wave profiles for the case kA ¼ 0.21 along the centerline at t ¼ 8T with and without the use of smoothing. The five-node smoothing technique was used every eight time steps for the case kA ¼ 0.21 and every four time steps for the case kA ¼ 0.31. It can be seen that the smoothing technique used in the present study retains good numerical stability. Fig. 9 shows the 3D profile at t ¼ 12T for the case kA ¼ 0.31, confirming that the 2D solution is accurately represented in the 3D tank.

v ¼ vð1Þ ¼

N X

bi eki z ki cos bi ;

(21)

i¼1

and the corresponding free surface elevation is given by Z ¼ Zð1Þ ¼

N X

ai cos bi ,

(22)

i¼1

where ai o i bi ¼ ; ki

bi ¼ ki ðx  x0 Þ  oi ðt  t0 Þ.

Here, ai and oi are the amplitude and frequency of the ith wave component, N the total number of wave components, and x0 and t0 correspond to the position and time where the predetermined largest wave crest occurs. For the non-linear focused wave group, the incident velocity is specified to be the second-order equivalent Stokes analytical velocity potential for deep water (Longuet-Higgins, 1963). Thus, v ¼ vð1Þ þ vð2Þ

(23)

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

5th-order theory

13T

10.5T

1.5

1

1

0.5

0.5 η/A

η/A

10T

1.5

0

5th-order theory

0

-0.5

-0.5

-1

-1

-1.5

12.5T

893

-1.5 -10 -5

0

5

10 15 x (m)

20

25

30

-10 -5

0

5

10 15 x (m)

20

25

30

Fig. 7. Free surface elevation spatial profiles for two steepnesses.

The corresponding second-order wave elevation is

1.5

Z ¼ Zð1Þ þ Zð2Þ

1

with

η/A

0.5

Z

0 -0.5 -1

-10

-5

0

5

10 x(m)

without smoothing

15

20

25

with smoothing

0.4

Z

0.2 0 -0.2 30 20 10

Y

0 -10

X

Fig. 9. 3D wave profile at t ¼ 12T as kA ¼ 0.31.

with ¼

¼

"

 1=2 pffiffiffiffi pffiffiffiffi2 ki  kj 2 ki kj  ðki kj Þ1=2 pffiffiffiffi pffiffiffiffi2 ki  kj  ki kj i¼1 j¼i ! pffiffiffiffi pffiffiffiffi2 2ðki kj Þ1=2 ki  kj þ ðki þ kj Þ cos bi cos bj þ pffiffiffiffi pffiffiffiffi2 ki  kj  ki kj ! #

N 1 X N X

ai aj

N 1 X i¼1

# 2ki kj ðoi  oj Þ bi bj ðoi  oj Þ2  gjki  kj j j¼iþ1 N X

"

 ejki kj jz ðki  kj Þ cosðbi  bj Þ.

sin bi sin bj

(26)

30

Fig. 8. Comparison of results with and without smoothing.

v

ð2Þ

ðki kj Þ1=2

-1.5

ð2Þ

(25)

(24)

In Eqs. (23) and (24), n(1) and Z(1) correspond to the linear terms in Eqs. (21) and (22), respectively. Firstly, we test the behaviour of the present numerical scheme by simulating a linear focused wave. For this simulation, the linear boundary value problem is solved at each time step. It is not necessary to provide a close match to a specific physical wave group, as we are here simply validating the model. In this example the focused wave is obtained by linear superposition of waves with periods T=2.8, 3.0, 3.2, 3.4, 3.6 and 3.8 s. It is assumed that the total wave amplitude A=Sai=0.2 m, and the wave components are taken to be of equal amplitude. It is convenient to define lmin as the shortest wavelength. The principal dimensions of the wave tank are then specified as length 6lmin, width 0.14lmin, and the depth of the wave source surface is truncated at 1.5lmin. The corresponding linear focal position and focal time are specified as x0=1.5lmin and t0=8Tmin, respectively. It is found that in this case the waves may be suitably dissipated by taking the lengths of the damping layers at each end of the tank as 1.5 times the minimum wavelength. The mesh density distributed along the tank and down the depth of the source surface is 16 elements per minimum wavelength. The time history of the simulated wave elevation at the focal position (x0, B/2), and comparison with the analytical solution, are shown in Fig. 10, where the wave elevation is normalized by A and the time is normalized by focal time t0. From the figure, it can be seen that the simulated results

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

894

agree well with those from the linear theory and the elevation reaches its largest value at t ¼ t0, with a symmetric distribution of wave profile about the focus time t0. The surface profile along the plane of symmetry down the tank at the focal time t ¼ t0 is shown in Fig. 11, where the x-coordinate is normalized by the focal distance x0.

1.5 Present Theory

1

η/A

0.5 0 -0.5 -1 -1.5 -0.4

-0.3

-0.2

-0.1

0 0.1 (t-t0) / t0

0.2

0.3

0.4

Fig. 10. Distribution of wave elevation at (x0, B/2).

1.5 Present Theory

1

0 -0.5 -1 -1.5 -1

-0.5

0 (x-x0) / x0

0.5

1

Fig. 11. Wave profile along the symmetry line at t ¼ t0.

P (0.33x0, B/2)

Q (1.67x0, B/2) 1.5

1.5 Present Theory

1

Present Theory

1 0.5 η/A

0.5 η/A

η/A

0.5

It can be seen that the corresponding elevation is largest at x ¼ x0 and the profile is symmetric about the focal position. Again we can see good agreement between the simulated results and linear theory. Fig. 12 shows the time histories of wave elevation at two points: P(0.33x0, B/2) and Q(1.67x0, B/2), respectively. The numerical results are in good agreement with the theory except for the leading transient waves, which are affected by the ramp function used for the numerical simulation. Moreover, these two time histories are antisymmetric with respect to each other, as the two points P and Q are symmetric about the focal position x0. Similar observations can be made concerning the behaviour of the wave profiles at times symmetric about the focal time. Fig. 13 illustrates details of the development of the wave group close to the focal time and focal location. Results are shown for the surface profile at t ¼ t0, t0+0.14 s and t00.14 s, and time history of the surface elevations at x ¼ x0, x1 ¼ x0+Dx and x2 ¼ x0Dx, respectively. Adjacent to the focus time and position, the crest elevations are reduced, and the surfaces are antisymmetric. These characteristics will be used to determine the location of the focal position for the fully non-linear cases investigated next. For non-linear focused wave simulation, a Stokes wave based on Eq. (23) and formed by the second-order interaction of the above waves is used as input at the source surface at every time step. The same conditions are specified for the computation as for the linear focused wave group except for the focal position and focal time. These are selected to be such that in the case of linear waves the focus positions and times would be x0 ¼ lmin, t0 ¼ 7Tmin. When a steep focused wave group evolves in the wave tank, however, there is a downstream shifting of the focal position relative to that occurring for small amplitude waves, with a corresponding increase in the time to focusing. This has been demonstrated both experimentally (Baldock et al., 1996; Johannessen and Swan, 2000) and theoretically (Longuet-Higgins, 1974; Kim et al., 1993) in finite water-depth. The non-linear shifts in the focal event depend on the amplitude A. The actual focal position x1 and focal time t1 will therefore differ from the linear values (x0 and t0), as will be shown in the following results.

0

0

-0.5

-0.5

-1

-1

-1.5

-1.5 -1

-0.5

0 (t-t0) / t0

0.5

1

-1

-0.5

0 (t-t0) / t0

Fig. 12. Distribution of wave elevations at two locations.

0.5

1

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

1.5

1.5 t0-0.14s

t0

x1

t0-0.14s

1

1

0.5

0.5 η/A

η/A

895

0

-0.5

-1

-1

-1

-0.5

0 (x-x0) / x0

0.5

1

x2

0

-0.5

-1.5

x0

-1.5 -0.25 -0.2 -0.15 -0.1 -0.05

0 0.05 (t-t0) / t0

0.1

0.15

0.2

0.25

Fig. 13. Wave focusing close to focal time and focal point.

Fig. 14 concerns the time histories of surface elevation at crest and trough focused positions x1 resulting from different input wave amplitudes. Four values of A are considered, defined as the linear sum of component wave amplitudes (A ¼ 0.2, 0.4, 0.6 and 0.8 m). In each figure, the non-linear surface elevation at the focal position is compared with the first-order analytical solution derived from Eq. (22) and the second-order analytical solution from Eq. (25). Since the linear and second-order analytical solutions are always focused at the linear focal time t0, they are here shifted to the actual focal time t1 of the non-linear simulations so that it is convenient to compare the peak wave amplitude of the analytical solutions with those of the numerical results. With the smallest input amplitude (Fig. 14(a)), the difference between the simulated results and linear solutions is small. As the amplitude of the individual wave components is very small, the non-linear wave–wave interactions are almost nonexistent. Fig. 14(b)–(d) show results for the larger input amplitudes. In these cases, the wave crest at the focus position becomes higher, while adjacent wave troughs become less deep. The opposite occurs for the wave trough at the focal position. One can also observe in Fig. 14(c) and (d) a significant difference before the focal time between the non-linear simulated results and the two analytical solutions for the larger input wave amplitude. This is related to the dependence of the celerities of the different wave components on wave amplitude, at orders of non-linearity above first- and second-order. The expected steeper wave envelopes for crest focused waves and broader wave envelopes for trough-focused waves are illustrated in the cases of increasing wave non-linearity. The results also confirm that the differences between focal times t1 and t0 become larger with the increase of input wave amplitude A. Non-linear features of the focused wave groups in Fig. 14 have also been investigated in the context of a spectral analysis. Fig. 15 compares the spectra of the numerically

simulated non-linear wave groups with those representing the linear and second-order analytical solutions at crest and trough focal points, respectively. (As before, these have been estimated by obtaining FFTs from the time histories of the wave elevations.) In all cases there is good agreement between numerical results and the linear and second-order analytical solutions for the spectrum in the range (0.26 Hzofo0.36 Hz). The figures also show that the non-linear effects from second-order difference and superharmonic waves, third-order super-harmonic waves—and even fourth-order super-harmonic waves in Fig. 15(d)— become increasingly significant as the input wave amplitude A is increased. For the crest and trough focused waves with the same input wave amplitude A, the corresponding spectra are seen to be very similar. Contributions to the wave elevation from linear and higher order terms are illustrated in Fig. 16 for the case A ¼ 0.4 m corresponding to Fig. 15(b). The time histories of crest and trough focal wave elevation from different frequency bands are obtained using the IFFT. The secondorder difference term corresponding to the long wave is shown in Fig. 16(a) and (b) plots the wave obtained from components in the linear range, showing good comparison between the numerical result and analytical solution. Fig. 16(c) shows the wave elevation from second-order superharmonic terms; and Fig. 16(d) illustrates the third-order super-harmonic terms The results confirm the obvious requirement that whereas at the focus time the first- and third-order contributions to the elevation are in anti-phase for the crest focused and trough focused waves, the secondorder contributions are in-phase: indeed the two plots in Fig. 16(c) are very similar. 5. Conclusions A time-domain numerical scheme has here been implemented to simulate linear and non-linear, regular and

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

896

1.5

1.5 1

crest

1 0.5 η /A

η /A

0.5 0

-0.5

-1

-1

-1.5 -0.25

0 (t-t0)/t0

-1.5 -0.25

0.25

crest

1 η /A

η /A

0

-1 0 (t-t0)/t0

-1.5 -0.25

0.25

0 t/t0

0.25

1.5

1.5 crest

1

trough

0.5 η /A

0.5 η/A

0.25

0

-1

-1.5 -0.25

0

0

-0.5

-0.5

-1

-1

-1.5 -0.25

0 (t-t0)/t0

-1.5 -0.25

0.25

1.5

1.5 crest

1

trough

0.5 η /A

0.5 0

0

-0.5

-0.5

-1

-1

-1.5 -0.25

0 t/t0

trough

-0.5

-0.5

1

0.25

0.5

0.5

1

0 (t-t0)/t0

1.5

1.5

η/A

0

-0.5

1

trough

0 (t-t0)/t0

0.25 Present

-1.5 -0.25

Linear theory

0 (t-t0)/t0 2nd-order theory

0.25

Fig. 14. Time histories of wave elevation at focus point x1 with different amplitude A (a) A ¼ 0.2 m, (b) A ¼ 0.4 m, (c) A ¼ 0.6 m and (d) A ¼ 0.8 m.

focused waves in a tank of infinite water-depth, by using a simple source higher order boundary element method. An image Green function is applied in the whole fluid domain so as to exclude the lateral boundaries, and transform the infinite domain to a finite domain. The linear and fully non-linear boundary value problems have

been solved separately, with appropriate boundary conditions. The resulting quadratic boundary integral equations have been solved within a time-stepping procedure. The instantaneous wave elevations and potentials on the free surface are updated using the fourth-order Runge–Kutta method. The incident waves are generated by means of a

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

70

70 crest

60

50 S(f)

S(f)

40 30

30 20

10

10 0 0

0.2

0.4

0.6 0.8 f(Hz)

1

1.2

1.4

0

0.2

0.4

0.6 0.8 f(Hz)

1

1.2

1.4

70

70 crest

60

50 S(f)

40 30

trough Present Linear 2nd-order

60

Present Linear 2nd-order

50 S(f)

40

20

0

40 30

20

20

10

10 0

0 0

0.2

0.4

0.6 0.8 f(Hz)

1

1.2

1.4

0

0.2

0.4

0.6 0.8 1 f(Hz)

1.2

1.4

70

70 crest

60

50 S(f)

40 30

trough Present Linear 2nd-order

60

Present Linear 2nd-order

50 S(f)

trough Present Linear 2nd-order

60

Present Linear 2nd-order

50

40 30

20

20

10

10

0

0 0

0.2

0.4

0.6 0.8 f(Hz)

1

1.2

1.4

0

70

0.2

0.4

0.6 0.8 f(Hz)

1

1.2

1.4

70 crest

60

50 S(f )

40 30

40 30

20

20

10

10

0

trough Present Linear 2nd-order

60

Present Linear 2nd-order

50 S(f)

897

0 0

0.2

0.4

0.6 0.8 f(Hz)

1

1.2

1.4

0

0.2

0.4

0.6 0.8 f(Hz)

1 1.2

1.4

Fig. 15. Comparisons of wave spectrum from fully non-linear results and analytical solutions by FFT (a) A ¼ 0.2 m, (b) A ¼ 0.4 m, (c) A ¼ 0.6 m and (d) A ¼ 0.8 m..

vertical plane of sources. The open boundaries are modelled by using damping layers. For the non-linear simulations, fifth-order Stokes theory for regular waves and second-order Stokes theory for focused wave groups

are used to define the incident waves which are fed through the source surface. Results from the linear and non-linear wave simulations have been compared with those calculated from the theory and good agreement obtained in

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

898

0.2

0.2 trough focused wave

crest focused wave

0.1 η/A

η/A

0.1 0

-0.1

-0.1 -0.2 -0.25

0 (t-t 0)/t 0

-0.2 -0.25

0.25

crest focused wave

η/A

η/A

0.5

0 -0.5

0 -0.5

-1

-1 0 (t-t 0)/t 0

-1.5 -0.25

0.25

0.2

0 (t-t 0)/t 0

trough focused wave

0.1 η/A

η/A

0.1 0 -0.1

0 -0.1

0 (t-t 0)/t 0

-0.2 -0.25

0.25

0.2

0 (t-t 0)/t 0

trough focused wave

0.1 η/A

0.1 η/A

0.25

0.2 crest focused wave

0

0 -0.1

-0.1 -0.2 -0.25

0.25

0.2 crest focused wave

-0.2 -0.25

0.25

trough focused wave

1

0.5

-1.5 -0.25

0 (t-t 0)/t 0

1.5

1.5 1

0

0 (t-t 0)/t 0 Present

-0.2 -0.25

0.25 Linear theory

0 (t-t 0)/t 0

0.25

2nd-order thoery

Fig. 16. Time histories of crest and trough focal wave elevation from different order wave components: time histories of wave elevation from (a) secondorder difference term, (b) linear term, (c) second-order super-harmonic term and (d) third-order superharmonic term.

cases of small amplitude. The physical characteristics of non-linear focused wave groups have been investigated using the resulting methodology. The present numerical scheme can be extended to the problems of fully non-linear wave–body interactions in very deep water.

Acknowledgements The authors gratefully acknowledge the financial support provided by the Program for Changjiang Scholars and Innovative Research Team in University (Grant no.

ARTICLE IN PRESS D.Z. Ning et al. / Ocean Engineering 35 (2008) 887–899

IRT0420), the Royal Society–NSFC China–UK Joint Project (Grant no. 15285) which facilitated exchanges between Dalian University of Technology and Oxford University and the National Natural Science Foundation of China (Grant No. 5063930, 50709005). References Baldock, T.E., Swan, C., Taylor, P.H., 1996. A laboratory study of nonlinear surface waves on water. Philosophical Transactions of the Royal Society of London 354A, 649–676. Beck, R.F., 1994. Time-domain computations for floating bodies. Applied Ocean Research 16, 267–282. Borthwick, A.G.L., Hunt, A.C., Feng, T., Taylor, P.H., Stansby, P.K., 2006. Flow kinematics of focused wave groups on a plane beach in the UK coastal research facility. Coastal Engineering 53 (12), 1033–1044. Brandini, C., Grilli, S.T., 2001. Modeling of freak wave generation in a 3D-NWT. In: Proceedings of the 11th Offshore and Polar Engineering Conference, vol III, Stavanger, Norway, pp. 124–131. Brebbia, C.A., Walker, S., 1980. Boundary Element Techniques in Engineering, Newnes-Butterworths. Brorsen, M., Larsen, J., 1987. Source generation of nonlinear gravity waves with the boundary integral equation method. Coastal Engineering 11, 93–113. Celebi, M.S., 2001. Nonliner transient wave–body interactions in steady uniform currents. Computer Methods in Applied Mechanics and Engineering 190, 5149–5172. Fenton, J.D., 1985. A fifth-order Stokes theory for steady waves. Journal of Waterway, Port, Coastal and Ocean Engineering, ASCE 111, 216–234. Ferrant, P., 1990. A coupled time and frequency approach for nonlinear wave radiation. In: Proceedings of the 18th Symposium on Naval Hydrodynamics, Michigan, USA, pp. 67–83. Ferrant, P., 1998. Fully nonlinear interactions of long-crested wave packets with a three-dimensional body. In: Proceeding of 22nd ONR Symposium on Naval Hydrodynamics, Washington, USA, pp. 403–415. Ferrant, P., Le Touze´, D., Pelletier, K., 2003. Nonlinear time domain models for irregular wave diffraction about offshore structures. International Journal of Numerical Methods in Fluids 43, 1257–1277. Fochesato, C., Grilli, S.T., Dias, F., 2007. Numerical modeling of extreme rogue waves generated by directional energy focusing. Wave Motion 44, 395–416. Hu, P., Wu, G.X., Ma, Q.W., 2002. Numerical simulation of nonlinear wave radiation by a moving vertical cylinder. Ocean Engineering 29 (14), 1733–1750. Johannessen, T.B., Swan, C., 2000. A laboratory study of the focusing of transient and directionally spread surface water waves. In: Proceedings of the Royal Society of London, vol. 457 A, pp. 971–1006. Kim, D.J., Kim, M.H., 1997. Wave–current–body interaction by a time-domain high-order boundary element method. In: Proceeding of the 7th International Offshore and Polar Engineering Conference (ISOPE-97), vol. 3, Honolulu, USA, pp. 107–115. Kim, C.H., Pierson, W.J., Tick, L.J., 1993. Extreme nonlinear wave loads on a vertical truncated circular cylinder in nonlinear irregular Stokeslike waves. In: Proceeding of the 3rd International Offshore and Polar Engineering Conference (ISOPE-93), vol. 3, Singapore, pp. 158–167. Kim, C.H., Cle´ement, A.H., Tanizawa, K., 1999. Recent research and development of numerical wave tank—a review. International Journal of Offshore and Polar Engineering 9 (4), 241–256. Lee, C.C., Liu, Y.H., Kim, C.H., 1994. Simulation of nonlinear waves and forces due to transient and steady motion of submerged sphere.

899

International Journal of Offshore and Polar Engineering 4 (3), 174–182. Li, H.B., Han, G.M., Mang, H.A., 1985. A new method for evaluating singular integrals in stress analysis of solids by the direct boundary element method. International Journal of Numerical Methods in Engineering 211, 2071–2075. Ligget, J.A., Liu, P.L., 1983. The Boundary Integral Equation Method for Porous Media Flow. George Allen & Unwin. Longuet-Higgins, M.S., 1963. The effect of nonlinearities on the statistical distributions in the theory of sea waves. Journal of Fluid Mechanics 17, 459–480. Longuet-Higgins, M.S., 1974. Breaking waves—in deep or shallow water. In: Proceedings of 10th ONR Symposium on Naval Hydrodynamics, Cambridge, MA, pp. 597–605. Longuet-Higgins, M.S., Cokelet, E.D., 1976. The deformation of steep surface waves on water. I. A numerical method of computation. In: Proceedings of the Royal Society of London, vol. 350 A, pp. 1–25. Newman, J.N., 1992. Approximation of free-surface Green functions. In: Martin, P.A., Whickham, G.R. (Eds.), Wave Asympotics. Cambridge University Press, pp. 107–142. Ning, D.Z., Teng, B., 2007. Numerical simulation of fully nonlinear irregular wave tank in three-dimension. International Journal for Numerical Methods in Fluids 53 (12), 1847–1862. Romate, J.E., 1989. The Numerical Simulation of Nonlinear Gravity Waves in Three Dimensions Using Higher-order Panel Method. Ph.D. Dissertation, University of Twente, Enschede, The Netherlands. Scorpio, S.M., Beck, R.F., Korsmeyer, F.T., 1996. Nonlinear water wave computations using a multipole accelerated desingularized method. In: Proceeding of 21st ONR Symposium on Naval Hydrodynamics, Trondheim, Norway, pp. 64–73. Taylor, P.H., Haagsma, I.J., 1994. Focusing of steep wave groups on deep water. In: Proceedings of the International Symposium: WavesPhysical and Numerical Modeling, Vancouver, Canada, pp. 862–870. Teng, B., Eatock Taylor, R., 1995. New higher-order boundary element methods for wave diffraction/radiation. Applied Ocean Research 17, 71–77. Teng, B., Ning, D.Z., 2006. Source generation of fully nonlinear waves in NWT with a 3D HOBEM. In: Proceedings of ISOPE PACOMS Conference, Dalian, China. Teng, B., Gou, Y., Ning, D.Z., 2006. A higher-order BEM for wavecurrent action on structures—direct computation of free term coefficient and CPV integrals. China Ocean Engineering 20 (3), 395–410. Tromans, P.S., Anaturk, A.R., Hagemeijer, A., 1991. A new model for the kinematics of large ocean waves-application as a design wave. In: Proceedings of 1st International Offshore and Polar Engineering Conference, vol. 3, Edinburgh, UK, pp. 64–71. Turnbull, M.S., Borthwick, A.G.L., Eatock Taylor, R., 2003. Numerical wave tank based on a s-transformed finite element inviscid flow solver. International Journal for Numerical Methods in Fluids 42, 641–663. Wang, C.Z., Khoo, B.C., 2005. Finite element analysis of two-dimensional nonlinear sloshing problems in random excitations. Ocean Engineering 32, 107–133. Wu, G.X., Ma, Q.W., Eatock Taylor, R., 1998. Numerical simulation of sloshing waves in a 3D tank based on a finite element method. Applied Ocean Research 20, 337–355. Xu, H., 1992. Numerical Study of Fully Nonlinear Water Waves in Three Dimensions. ScD Dissertation, Department of Ocean Engineering, Massachusetts Institute of Technology, Cambridge, USA. Xu, Q., Pawlowski, J.S., Baddour, R.E., 1993. Simulation of nonlinear irregular waves by Green–Naghdi theory. In: Proceedings of 8th International Workshop on Water Waves and Floating Bodies, Newfoundland, Canada, pp. 165–167. Zienkiewicz, O.C., Morgan, K., 1983. Finite Elements and Approximation. Wiley.