Michael Bestehorn¹, Peder A. Tyvand², Thomas Michelitsch^3
1Department of Physics, Brandenburg University of Technology, Cottbus, Germany.
²Department of Mathematical Sciences and Technology, Norwegian University of Life Sciences, Ås, Norway.
³CNRS, Institut Jean le Rond d’Alembert, Sorbonne Université, Paris, France.
DOI: 10.4236/jamp.2019.71007   PDF    HTML   XML   649 Downloads   1,310 Views   Citations

Abstract

Starting from the 2D Euler equations for an incompressible potential flow, a dimension-reduced model describing deep-water surface waves is derived. Similar to the Shallow-Water case, the z-dependence of the dependent variables is found explicitly from the Laplace equation and a set of two one- dimensional equations in x for the surface velocity and the surface elevation remains. The model is nonlocal and can be formulated in conservative form, describing waves over an infinitely deep layer. Finally, numerical solutions are presented for several initial conditions. The side-band instability of Stokes waves and stable envelope solitons are obtained in agreement with other work. The conservation of the total energy is checked.

Keywords

Hydrodynamics, Ocean Waves, Deep-Water Waves, Numerical Solutions, Fractal Derivatives

Share and Cite:

1. Introduction

The procedure for deriving reduced equations by expressing the dependence of the dependent variables on a certain, geometrically distinguished spatial coordinate analytically is widespread in nonlinear pattern formation and wave dynamics. As a first attempt, one can consider the Shallow-Water or Saint-Venant equations [1] , which are systematically derived from the hydrodynamic basic equations for an inviscid potential flow by solving iteratively the Laplace equation for the velocity potential function. Thereby, all expansions are controlled by a geometry parameter

$\delta =\frac{d}{\mathcal{l}}\ll 1$ (1)

where d is the layer’s depth and $\mathcal{l}$ a typical lateral length scale (e.g. wave length). A similar expansion can be done for highly viscous problems and the so-called Thin-Film equation is obtained which describes the dynamical behavior of the shape of the film’s surface [2] [3] [4] . For convection problems, a projection onto certain z-dependent modes allows for the reduction to 2D systems where only the horizontal coordinates occur. Standard equations in normalized form like the Allen-Cahn, the (complex) Ginzburg-Landau or the Swift-Hohenberg equations can be derived, partially by heuristic arguments concerning at least their nonlinear parts [5] .

Stokes [6] initiated the theory of nonlinear deep-water waves. Three branches of development for strongly nonlinear theory have grown from seminal papers. 1) The highest Stokes wave [7] . 2) The instability of Stokes waves [8] . 3) Breaking of deep-water waves [9] .

The early nonlinear theories in 2D have now grown into full 3D theories. However, a long-standing problem is that of dimensional reduction of the nonlinear 3D problem to a 2D analysis, as it can be done straightforwardly for the case of shallow water. For deep water waves, the Hamiltonian formulation [10] [11] provides a nonlinear system for two two-dimensional field variables, namely surface elevation and potential at the surface. To compute the velocity normal to the surface, the solution of the underlying 3D-Laplace equation is achieved by an expansion with respect to powers of the surface elevation into spectral function. These algorithms have been named ‘higher order spectral methods’ (HOSM) and go back to the work of West et al. [12] and Dommermuth and Yue [13] . Another attempt is the systematic derivation of an equation describing the slowly varying envelope of a given monochromatic wave train. This equation has the form of a nonlinear cubic Schrödinger equation [10] [14] and reflects the side-band instability of Stokes waves as well as envelope solitons as localized stable solutions.

Contrary to the Shallow-Water or Thin-Film equations, no intrinsic length scale exists for the deep-water case. The only available scales are given by the solutions sought for, namely the wave length $\lambda$ (if the waves are nearly coherent) and the amplitude A of those waves. In deep-water waves, the ratio of amplitude to wave length is normally rather small. Thus, using the wave steepness

$s=Ak$ (2)

with $k=2\text{π}/\lambda$ for expanding the nonlinearities seems to be quite natural.

The present paper offers a consistent scheme of reducing the spatial dimension of nonlinear deep-water wave theory by introducing the concept of fractional derivatives along the free surface. Starting from the 2D Euler equations for an incompressible potential flow, a one-dimensional model describing deep-water surface waves is derived. Similar to the Shallow-Water case, the z-dependence of the dependent variables is determined from the Laplace equation and a set of two nonlinear equations for the surface velocity and the surface elevation remains. The system is nonlocal and can be formulated in conservative form, describing waves over an infinitely deep layer. It contains nonlinearities of the lowest (quadratic) order in velocity and surface elevation but can be systematically extended to higher orders. Numerical solutions are presented for different initial conditions. Coherent wave trains prove to be unstable due to the Benjamin-Feir instability [8] and localized solutions in form of envelope solitons are obtained. Finally, a spatially localized and temporally oscillating surface pressure variation is used to act as a wave-maker. The propagation of different wave ensembles is studied and linear and nonlinear results are compared to show the importance of the nonlinear terms.

The conservation of the total energy can be considered as a quality test for both the derived system and the numerical method applied for its solutions. We show that the energy is conserved within a few percents which give us confidence in our derivation and results.

The model can be extended to the case of a 3D potential flow, resulting in a set of three two-dimensional evolution equations.

2. The Model

2.1. Potential Flow

For the sake of simplicity we restrict the treatment first to two dimensions, $v=\left(v_x\mathrm{,}{v}_z\right)$ . We are looking for solutions of the non-dimensional Euler eqs. for an incompressible fluid in the half space $-\infty <z\le h\left(x,t\right)$ , see Figure 1:

${\partial }_t{v}_x+v_x{\partial }_x{v}_x+v_z{\partial }_z{v}_x=-{\partial }_{x}P$ (3)

${\partial }_t{v}_z+v_x{\partial }_x{v}_z+v_z{\partial }_z{v}_z=-{\partial }_{z}P-1$ (4)

${\partial }_x{v}_x+{\partial }_z{v}_z=0$ (5)

where time is scaled by

$\tau =\sqrt{\frac{\mathcal{l}}{g}}$ (6)

and lengths with an arbitrary $\mathcal{l}$ which will be identified later (sect. 3.2) with the characteristic wave length of the initial condition. We assume that the flow field can be derived from a potential

$v=\nabla \Phi$ (7)

which then, due to (5) must be a solution of the Laplace eq.

${\nabla }^2\Phi =0$ (8)

which also implies

${\nabla }^{2}v=0.$ (9)

Let

$u_h\left(x,t\right)=v_x\left(x,z=h\left(x,t\right),t\right)$ (10)

be the horizontal velocity component at the free surface located at $z=h\left(x,t\right)$ .

From (3) we derive

${\partial }_t{u}_h=-\frac{1}{2}{\partial }_x{\left(u_h\right)}^2-{{\partial }_{x}P|}_{z=h}.$ (11)

Equation (11) serves as a boundary condition (b.c.) for (9) at $z=h\left(x,t\right)$ . The surface function h is determined by the kinematic b.c.

${\partial }_{t}h={v_z|}_{z=h}-u_h{\partial }_{x}h$ (12)

For the following treatment, it is of advantage to formulate (12) also in conservative form

${\partial }_{t}h=-{\partial }_{x}S$ (13)

with the flux

$S\left(x\mathrm{,}t\right)\mathrm{<mo>\; =\; </mo>}{\displaystyle {\int }_{-\infty }^{h\left(x\mathrm{,}t\right)}}\text{d}z{v}_x\left(x\mathrm{,}z\mathrm{,}t\right)\mathrm{.}$ (14)

For an infinitely deep layer, the asymptotic b.c.

$v=0\text{for}z\to -\infty$ (15)

must hold. Hence the general solution of (9) reads

$v_x={\displaystyle \underset{k}{\sum }}{u}_k\left(t\right){\text{e}}^{\left|k\right|z}{\text{e}}^{ikx}+c.c.$ (16)

with c.c. as the complex conjugate.

2.2. Nonlocal Expansion

Considering (11) and (13) as evolution Equations, the original 3D (here 2D) problem is reduced to a 2D (here 1D) system in the horizontal coordinates x,y (here x) only. However, to evaluate the flux (14) one needs to know $v_x\left(x\mathrm{,}z\mathrm{,}t\right)$ . Taking (16), Equation (14) reads

$S\left(x,t\right)={\displaystyle \underset{k}{\sum }}\frac{{\text{e}}^{\left|k\right|h\left(x\mathrm{,}t\right)}}{\left|k\right|}{u}_k\left(t\right){\text{e}}^{ikx}+c.c.$ (17)

To determine the amplitudes $u_k$ from $u_h$ , we evaluate (16) at $z=h$ and expand the exponential function up to a given order $h^N$ :

$u_h={\displaystyle \underset{k}{\sum }}{u}_k\left(t\right){\displaystyle \underset{n=0}{\overset{N}{\sum }}}\frac{{\left(\left|k\right|h\right)}^n}{n!}{\text{e}}^{ikx}+c.c.$ (18)

Next we introduce the fractional differential operator of order $h^N$ :

${\hat{L}}_N\left(h\right)={\displaystyle \underset{n=0}{\overset{N}{\sum }}}\frac{h^n{\left(-{\partial }_{xx}\right)}^{n/2}}{n!}$ (19)

with the nonlocal operator ${\left(-{\partial }_{xx}\right)}^{n/2}$ defined as

${\left(-{\partial }_{xx}\right)}^{n/2}f\left(x\right)={\displaystyle {\int }_{-\infty }^{\infty }}{\left|k\right|}^n{\tilde{f}}_k{\text{e}}^{ikx}\text{d}k$ (20)

for every function $f\left(x\right)$ via its Fourier transform ${\tilde{f}}_k$ . We note that (20) is one of the definitions of the fractal Laplacian, here in 1D [15] . Now Equation (18) can be written as

$u_h={\hat{L}}_N\left(h\right){\displaystyle \underset{k}{\sum }}{u}_k\left(t\right){\text{e}}^{ikx}+c.c.$ (21)

(in the 3D case, ${\partial }_{xx}$ must be replaced by the horizontal Laplacian ${\partial }_{xx}+{\partial }_{yy}$ ).

Using the same technique, S from (17) can be expressed as

$S_N={\hat{L}}_N\left(h\right){\left(-{\partial }_{xx}\right)}^{-1/2}{\displaystyle \underset{k}{\sum }}{u}_k\left(t\right){\text{e}}^{ikx}+c.c.$ (22)

leading finally to

$S_N\left(x,t\right)={\hat{L}}_N\left(h\right){\left(-{\partial }_{xx}\right)}^{-1/2}{\hat{L}}_N^{-1}\left(h\right)u_h\left(x,t\right)$ (23)

where ${\hat{L}}_N^{-1}$ is the operator inverse to ${\hat{L}}_N$ and is computed from

${\hat{L}}_n^{-1}{\hat{L}}_n=1+O\left(h^{n+1}\right).$ (24)

Thus, the equations (11) and (13) are closed by the relation (23). Note that ${\left(-{\partial }_{xx}\right)}^{n/2}$ is a nonlocal operator only for odd n. For n even it can be expressed by spatial derivatives of order n in position space. Moreover, assuming periodic lateral boundary conditions, it is self-adjoint for all n.

2.3. Second Order Model

In the following we shall restrict ourselves to $N=1$ , resulting in a model of second order in the dependent variables $h\left(x\mathrm{,}t\right)\mathrm{,}{u}_h\left(x\mathrm{,}t\right)$ . With

${\hat{L}}_1=1+h\hat{D},$ (25)

its inverse

${\hat{L}}_1^{-1}=1-h\hat{D},$ (26)

and the self-adjoint operator

$\hat{D}\equiv {\left(-{\partial }_{xx}\right)}^{1/2}\mathrm{,}$ (27)

Equation (23) turns into

$S_1={\hat{D}}^{-1}{u}_h-{\hat{D}}^{-1}\left(h\hat{D}{u}_h\right)+h{u}_h$ (28)

where only terms up to the second order in $u_h\mathrm{,}h$ are included.

In (11), the pressure gradient evaluated at $z=h$ turns into

${{\partial }_{x}P|}_{z=h}={\partial }_x{P}_h-{{\partial }_{z}P|}_{z=h}{\partial }_{x}h={-{\partial }_{z}P|}_{z=h}{\partial }_{x}h$ (29)

where $P_h$ is the pressure along the free surface and assumed to be constant (external pressure). To find ${\partial }_{z}P$ we evaluate (4) at the surface and include only linear terms in $v$ , leading to

${{\partial }_{z}P|}_{z=h}=-1-{\partial }_t{w}_h$ (30)

with the vertical velocity component

$w_h\left(x,t\right)=v_z\left(x,z=h\left(x,t\right),t\right)$ (31)

along the surface. $w_h$ is related to $u_h$ by (5) according to

$w_h=-{\partial }_x{\hat{D}}^{-1}{u}_h$ (32)

where again only linear terms are considered. Finally, we have

${{\partial }_{x}P|}_{z=h}={\partial }_{x}h\left(1-{\partial }_x{\hat{D}}^{-1}{\partial }_t{u}_h\right),$ (33)

where the first term on the right-hand side corresponds to the hydrostatic pressure, the second one to the dynamic part. Using (33), Equation (11) can be cast into the form

$\left[1-\left({\partial }_{x}h\right){\partial }_x{\hat{D}}^{-1}\right]{\partial }_t{u}_h=-{\partial }_x\left[\frac{1}{2}{u}_h^2+h\right].$ (34)

2.3.1. Model

After inversion of the operator on the left-hand side of (34), the complete system up to the second order reads

${\partial }_t{u}_h=-{\partial }_x\left[\frac{1}{2}{u}_h^2+h\right]+\left({\partial }_{x}h\right)\hat{D}h$ (35)

${\partial }_{t}h=-{\partial }_x{S}_1$ (36)

with $S_1$ from (28).

2.3.2. Linear Waves

Taking only linear terms into account one finds from ((35), (36))

${\partial }_t{u}_h=-{\partial }_{x}h$

${\partial }_{t}h=-{\partial }_x{\hat{D}}^{-1}{u}_h$ (37)

or, eliminating $u_h$ , the linear wave equation

${\partial }_{tt}h={\partial }_x{\hat{D}}^{-1}{\partial }_{x}h=-\hat{D}h$ (38)

where we have used the identity ${\partial }_{xx}=-{\left(\hat{D}\right)}^2$ . Equation (38) possesses the well-known short-wave dispersion relation

$\omega ={\left|k\right|}^{1/2}.$ (39)

2.3.3. Two Horizontal Dimensions

If we add the second horizontal dimension we are able to describe a 3D flow which, however, is still potential. Extending the surface velocity to a 2D vector

$u_h\left(x,y,t\right)=\left({\left(u_h\right)}_x,{\left(u_h\right)}_y\right)$ (40)

the derivation of the reduced system is straightforward and it reads

${\partial }_t{u}_h=-{\nabla }_2{S}_u+\left({\nabla }_{2}h\right)\hat{D}h$ (41)

${\partial }_{t}h=-{\nabla }_2\cdot S_1$ (42)

$S_u=\frac{1}{2}{\left(u_h\right)}^2+h$ (43)

$S_1={\hat{D}}^{-1}{u}_h-{\hat{D}}^{-1}\left(h\hat{D}{u}_h\right)+h{u}_h$ (44)

with the 2D nabla operator ${\nabla }_2=\left({\partial }_x,{\partial }_y\right)$ .

2.4. Similarities and Differences with HOSM

Contrary to the higher-order spectral methods (HOSM) presented in [12] [13] , we use the surface elevation and the surface velocity instead of the surface potential. With the kinematic b.c. (36) formulated by the help of the flux S, this results in model equations that can be written in conservative form. Due to the self-adjointness of $\hat{D}$ it can be easily seen that ${\displaystyle {\int }_{\Gamma }}\text{d}x{u}_h\mathrm{,}{\displaystyle {\int }_{\Gamma }}\text{d}xh$ are conserved in time, also for the numerical scheme, if the derivatives with respect to x are expressed by simple central differences (see sect. 3.1 below). On the other hand, our method can be formulated in a closed form in position space and is, at least up to second order, faster and easier to implement numerically, allowing for larger aspect ratios and an effective numerical realization in two horizontal dimensions. Eventually, Equations ((35), (36)) can be considered as a systematic expansion with respect to the wave steepness $s=kh$ . With the self-consistent assumption $u_h=O\left(s\right)$ , our model includes all terms up to the order $s^2$ . Hence it should deliver accurate results up to second-order stokes waves, which corresponds to $s_m\approx 0.16$ [16] , about 36 percents of the steepest possible Stokes wave with $s=0.44$ . Then the maximal wave amplitude of a wave with $k=2\text{π}$ (wave length $\lambda =1$ ) should not be much higher than $h_{\max }\approx 0.026$ .

3. Numerical Method and Results

3.1. Numerical Method

We concentrate on the 2^nd order model ((35), (36)) with x in the domain $\left[\mathrm{0,}\Gamma \right]$ and periodic lateral b. c. $u_h\left(0\right)=u_h\left(\Gamma \right),h\left(0\right)=h\left(\Gamma \right)$ . To evaluate (28), the expressions

$\hat{D}f\left(x\right)\mathrm{,}{\hat{D}}^{-1}f\left(x\right)$ (45)

must be computed. The operator introduced in (20) thus reads ${\hat{D}}^n={\left(-{\partial }_{xx}\right)}^{n/2}$ . It is nonlocal for odd n and (45) can be written as

${\hat{D}}^{n}f\left(x\right)={\displaystyle {\int }_0^{\Gamma }}\text{d}{x}^{\prime }{G}^{\left(n\right)}\left(x-x^{\prime }\right)f\left(x^{\prime }\right)$ (46)

with the Green’s function

$G^{\left(n\right)}\left(x\right)=\frac{1}{\Gamma }{\displaystyle \underset{\mathcal{l}=-\infty }{\overset{\infty }{\sum }}}{\left|k_{\mathcal{l}}\right|}^n{\text{e}}^{i{k}_{\mathcal{l}}x}$ (47)

and $k_{\mathcal{l}}=2\text{π}\mathcal{l}/\Gamma$ . Taking a finite difference scheme for the numerical integration of the system, the variable x is discretized as

$x_j=j\Delta x,\Delta x=\Gamma /N,j=0,\cdots ,N-1$

where N is the number of mesh points and $f_N=f_0$ with $f_i\equiv f\left(x_i\right)$ . The discrete form of (46) thus reads

${\hat{D}}^n{f}_j={\displaystyle \underset{\mathcal{l}=0}{\overset{N-1}{\sum }}}{G}_{j\mathcal{l}}^{\left(n\right)}{f}_{\mathcal{l}}$ (48)

with

$G_{j\mathcal{l}}^{\left(n\right)}=\frac{1}{N}{\displaystyle \underset{m=-N/2}{\overset{N/2}{\sum }}}{\left|k_m\right|}^n{\text{e}}^{i{k}_m\left(x_j-x_{\mathcal{l}}\right)}$ (49)

which would result in time-consuming computations of convolutions when evaluating (28). Hence it is much more effective to compute (48) in Fourier space where it simply reads

${\hat{D}}^n{\tilde{f}}_{\mathcal{l}}={\left|k_{\mathcal{l}}\right|}^n{\tilde{f}}_{\mathcal{l}}$ (50)

where ${\tilde{f}}_{\mathcal{l}}$ is the discrete Fourier transform of $f_i$ that can be obtained numerically by any standard Fast-Fourier transform (FFT) [17] .

To evaluate S according to (28) it is therefore necessary to compute $\hat{D}{u}_h$ and ${\hat{D}}^{-1}{u}_h$ in Fourier space, then perform the products $h\hat{D}{u}_h$ and $h{u}_h$ in real space and finally compute ${\hat{D}}^{-1}h\hat{D}{u}_h$ in Fourier space again. In addition for (35), $\hat{D}h$ has to be determined in Fourier space. Thus, seven FFTs (three forward and four backward) have to be performed at each time step.

The derivatives with respect to x occurring in ((35), (36)) are approximated with centered differences

${{\partial }_{x}f|}_{x_i}\to \frac{f_{i+1}-f_{i-1}}{2\Delta x}$ (51)

what ensures the conservation of ${\displaystyle {\sum }_i}{f}_i$ .

Finally, time integration is achieved by an explicit 4^th-order Runge-Kutta algorithm [18] .

3.2. Results

3.2.1. Benjamin-Feir Instability of Slightly Disturbed Coherent Wave Trains

The initial conditions are based on surface elevations in form of a traveling wave with given amplitude and wave number $k_0$ . If for the spatial scaling introduced in (6) the wave length is chosen, we can put $k_0=2\text{π}$ . The x-dimension is discretized with $N=2048$ mesh points, the step sizes used are $\Delta x=60/N\approx 0.029$ , $\Delta t=0.5\times {10}^{-3}$ , resulting in a side length of $\Gamma =60$ and allowing for 60 initial wave lengths. Each wave length is resolved with about 33 mesh points.

As first shown in [8] , coherent Stokes waves with a single frequency and wave length are unstable with respect to the side-band instability. The time series in Figure 2 presents the evolution of a 2^nd-order Stokes wave [19] traveling to the right side as initial condition:

$h\left(x,0\right)=A\cos \left(k_{0}x\right)+\frac{1}{2}{k}_0{A}^2\left(1+\cos \left(2k_{0}x\right)\right)+\epsilon f(x)$

$u_h\left(x,0\right)=\sqrt{k_0}\left(A\cos \left(k_{0}x\right)+\frac{3}{4}{k}_0{A}^2\left(1+\cos \left(2k_{0}x\right)\right)+\epsilon f\left(x\right)\right)$ (52)

with

$A=0.0075,k_0=2\text{π},{\omega }_0=\sqrt{k_0}.$

To trigger the side band (Benjamin-Feir, BF) instability [20] , weak disturbances of the form

$f\left(x\right)=\cos \left(\left(k_0+\delta k\right)x\right)+\cos \left(\left(k_0-\delta k\right)x\right)$ (53)

with $\epsilon =0.01A$ and $\delta k=2\text{π}/\Gamma$ have been added. The speed and wave length of the (undisturbed) traveling wave follows as

$c_0={\omega }_0/k_0+O\left(A^2\right)=k_0^{-1/2}+O\left(A^2\right)\approx 0.4,$ (54)

what yields a Courant number of

$c_0\Delta t/\Delta x\approx \mathrm{0.007.}$

The disturbances (39) belong to neighbored waves of $k_0$ in Fourier space. The absolute values of the Fourier amplitudes over time for the carrier wave ( $k_0$ ) as well as for its five neighbors ( $k_0\pm \delta k\mathrm{,}{k}_0\pm 2\delta k$ ) are shown in Figure 3. In the beginning, the side bands grow exponentially, taking the energy from the carrier wave that decreases.

3.2.2. Envelope Solitons

If a coherent wave is enveloped by a slowly varying function $\Phi \left(x\mathrm{,}t\right)$ according to

$h\left(x,t\right)=\mathrm{Re}\left[\Phi \left(x,t\right){\text{e}}^{i\left(k_{0}x+\omega t\right)}\right],$ (55)

Zakharov [10] showed that for small but finite amplitudes an equation for the complex valued $\Phi \left(x\mathrm{,}t\right)$ can be derived, having the form of the nonlinear Schrödinger equation (NLS):

$i{\partial }_t\Phi =\frac{1}{2}{\omega }^{″}{\partial }_{\xi \xi }\Phi +\nu {\left|\Phi \right|}^2\Phi$ (56)

with $\xi =x-{\omega }^{\prime }t$ and

${\omega }^{\prime }={\text{d}\omega /\text{d}k|}_{k_0}=\frac{1}{2}{k}_0^{-1/2},{\omega }^{″}={{\text{d}}^2\omega /\text{d}{k}^2|}_{k_0}=-\frac{1}{4}{k}_0^{-3/2}.$ (57)

The nonlinear coefficient $\nu$ depends on the ratio of the water depth to the wave length of the carrier wave [11] and reads for infinite depth

$\nu =-\frac{1}{2}{k}_0^{5/2}.$ (58)

Equation (56) has the well-known localized hyperbolic secant solutions

$\Phi \left(\xi \mathrm{,}t\right)=\frac{A}{\cosh \left(\kappa \xi \right)}{\text{e}}^{i\Omega t}$ (59)

with the relation between reciprocal width and amplitude

$\kappa =\sqrt{\frac{\nu }{{\omega }^{″}}}A=\sqrt{2}{k}_0^{2}A$ (60)

and $\Omega =-{\kappa }^2{\omega }^{″}/2=k_0^{5/2}{A}^2/4$ .

To see how envelope solitons behave in our model, we take as initial condition

$h\left(x,0\right)=\frac{A}{\cosh \left(\kappa \left(x-x_0\right)\right)}\sin \left(k_{0}x\right),u_h\left(x,0\right)={\omega }_{0}h\left(x,0\right).$ (61)

Figure 4 shows snapshots from the evolution for

$A=0.006,\kappa =1/3,x_0=0.2\Gamma ,k_0=2\text{π},{\omega }_0=\sqrt{k_0}$

with $N=2048$ and $\Gamma =60$ . The envelope soliton travels thereby with the group speed ${\omega }^{\prime }$ as given in (57) which is half of the phase speed ${\omega }_0/k_0$ of the carrier wave.

To confirm the agreement of our model and its numerical solutions with the results of the weakly nonlinear NLS, we also studied an initial condition as (61), but with an amplitude and width of the envelope that does not satisfy the relation (60). Taking the same $\kappa =1/3$ but $A=0.003$ , the package flows apart and finally extends over the whole domain, Figure 5.

Finally we show in Figure 6 the interaction of two colliding envelope solitons, generated by the initial condition

$h\left(x,0\right)=\left[\frac{A_1}{\cosh \left({\kappa }_1\left(x-x_1\right)\right)}+\frac{A_2}{\cosh \left({\kappa }_2\left(x-x_2\right)\right)}\right]\sin \left(k_{0}x\right)$ (62)

$u_h\left(x,0\right)=\left[\frac{A_1{\omega }_0}{\cosh \left({\kappa }_1\left(x-x_1\right)\right)}-\frac{A_2{\omega }_0}{\cosh \left({\kappa }_2\left(x-x_2\right)\right)}\right]\sin \left(k_{0}x\right)$

with

$A_1=0.006,A_2=0.009,{\kappa }_1=1/3,{\kappa }_2=1/2,x_1=0.2\Gamma ,x_2=0.8\Gamma .$

The solitons travel in opposite direction and interact without changing their shapes strongly. Due to the periodic b.c. these interactions repeat after equidistant time intervals of $T=\Gamma /2{\omega }^{\prime }\approx 150$ .

3.2.3. Conservation of Energy

It is interesting to check the temporal variation of the energy during the runs. For an inviscid fluid, the total energy

$E=E_p+E_k$ (63)

with the potential energy

$E_p=\frac{1}{2}{\displaystyle {\int }_0^{\Gamma }}\text{d}x{h}^2$ (64)

and the kinetic energy

$E_k=\frac{1}{2}{\displaystyle {\int }_0^{\Gamma }}\text{d}x{\displaystyle {\int }_0^h}\text{d}z\left[v_x^2+v_z^2\right]$ (65)

must be conserved. For the latter one, one finds up to quadratic order the expression

$E_k=\frac{1}{2}{\displaystyle {\int }_0^{\Gamma }}\text{d}x{u}_h{\hat{D}}^{-1}{u}_h\mathrm{.}$ (66)

The conservation of (63) can serve as a quality feature for the numerical solutions as well as for the accuracy of the derived system itself.

Figure 7, Figure 8 show the three energies (63)-(65) over time for the series depicted in Figure 2 and Figure 6, respectively. In Figure 8, potential as well as kinetic energy varies strongly during the collision phases but the sum (63) remains fairly constant. Table 1 shows the standard deviations normalized with the square mean defined as

$\sigma \left(f\right)=\sqrt{\frac{\langle f{\left(t\right)}^2\rangle -{\langle f\left(t\right)\rangle }^2}{\langle f{\left(t\right)}^2\rangle }}$ (67)

with the mean

$\langle \mathrm{...}\rangle =\frac{1}{T}{\displaystyle {\int }_0^T}\text{d}t\mathrm{...,}$ (68)

$T=7000$ and f as one of the three energy functions.

Note that all computations are only up to second order which then applies also for the energy conservation so that

$d_{t}E=O\left({\xi }^3\right)$ (69)

where $\xi$ stands for h, $u_h$ and combinations.

3.2.4. High Waves Caused by Dispersive Focusing of an Unstable Stokes Wave

Waves on deep water are generated by wind and have a broad spectrum of wave lengths. Since their phase and group speed is proportional to ${\omega }^{-1}$ , longer waves are faster than shorter ones and, if they propagate in the same direction, may catch up and enhance the total amplitude. This is already a linear effect but can be amplified by nonlinearities what then leads to even higher wave heights $H\approx 2h$ . In this way, very large heights can be obtained, which is one of the mechanisms responsible for the occurrence of freak or rogue waves [14] and named “dispersive focusing” or “dispersive enhancement”.

To classify waves with respect to their height, the “abnormality index” AI

$\text{AI}=H_m/H_s$ (70)

is introduced, where $H_m$ is the height of the wave and $H_s$ the significant wave height. The significant wave height is defined as the average height of the upper third of the wave height distribution (over time or space). Another way to compute $H_s$ is to use the root of the zeroth order moment of the spectrum

$H_{m0}=4{\left[\frac{1}{\Gamma }{\displaystyle {\int }_0^{\Gamma }}\text{d}x{h}^2\left(x\right)\right]}^{1/2}\mathrm{.}$ (71)

Then, a freak wave can be defined by AI > 2. The probability for freak waves depends strongly on the initial wave pattern and on the dimensionality of the system. In 2D, only two propagation directions are possible and freak waves should occur quite often. Figure 9 shows the evolution of an unstable Stokes wave as (52), but for a larger aspect ratio and a higher grid resolution $N=2^{14}=16384$ at $t=14200$ . A freak wave shows up with AI ≈ 2.8. In Figure 10 we plot AI over time, the circle marks the situation shown in Figure 9.

Figure 11 shows the the absolute values of the Fourier amplitudes of $h\left(x\right)$ depicted in Figure 9. A scaling law according to

$\left|h_k\right|\sim \exp \left(-\alpha \left|k\right|\right)$ (72)

with $\alpha \approx 2/15$ can be seen over wide regions, in agreement with the findings of Zakharov et al. [21] .

3.2.5. Waves Generated by Localized Pressure Forcing

In this section we consider an initially quiescent flat layer with $h=u_h=0$ at

$t=0$ on which waves are generated by a spatially localized surface pressure of the form

$P_h\left(x\mathrm{,}t\right)=P_0\left(t\right)\exp \left(-x^2/{\gamma }^2\right)\mathrm{.}$ (73)

Inclusion of (73) modifies (35) to

${\partial }_t{u}_h=-{\partial }_x\left[\frac{1}{2}{u}_h^2+h+P_h\right]+\left({\partial }_{x}h\right)\hat{D}h.$ (74)

If the amplitude $P_0$ varies harmonically, $P_0\left(t\right)=A\sin \left(\Omega t\right)$ (73) acts as a wave maker and produces waves with wave length $\lambda =2\text{π}/k={\text{2π}/\Omega }^2$ , traveling symmetrically to both sides. To simulate the situation where waves with shorter wave lengths are followed by longer (faster) waves, we perform a run with $\gamma =20\Delta x$ and $\Delta x=\Gamma /N=0.0146$ . Frequency and amplitude of the pressure are decreasing linearly in time, according to

$\Omega \left(t\right)=\max \left({\omega }_0-\frac{1}{4}\frac{t}{x_s},0\right),A\left(t\right)=\frac{\Omega \left(t\right)}{{\omega }_0}{A}_0$ (75)

with $A_0=0.005,{\omega }_0=\sqrt{\text{2π}},x_s=120$ . Two runs have been performed, one with the nonlinear model ((74), (36)), the other with the same system but only keeping linear terms in $h\mathrm{,}{u}_h$ . From Figure 12 it is evident that the nonlinearities amplify the amplitudes and provide significantly higher waves. Maximal focusing is obtained after $t\approx 590$ at $x_s$ , see Figure 13.

4. Conclusions

We presented the derivation and numerical solutions of dimension-reduced deep-water equations. The system can be considered as a systematic expansion

of the basic hydrodynamic equations with respect to powers of the wave steepness. We discussed the model in the lowest nonlinear order and showed numerical 1D-solutions where the Benjamin-Feir instability of coherent waves turned out and the behavior of envelope solitons has been studied. All results are in good agreement with previous work. The envelope solitons coincide well with solutions of the nonlinear Schrödinger equation both in amplitude and width. Our model is not based on a certain carrier wave and no spatial direction along the surface is distinguished. The formulation is rotationally invariant in the horizontal dimension(s) and allows therefore also for the study of colliding wave trains.

From the numerical side, the great advantage of the present model is its relatively straightforward and effective numerical realization. Due to the dimension reduction, CPU-times of the runs shown are in the range of hours on normal PCs or Laptops. Our system can be extended both to higher order nonlinearities as well as to three-dimensional potential flows, resulting in a reduced 2D set. Thus, the numerical examination of 2D envelope solitons occurring on the surface of a 3D potential flow should become feasible. This work is currently in progress and will be presented in a following paper.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Saint-Venant, A.J.C. (1871) Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l’introduction de marées dans leurs lits. Comptes Rendus des Séances de Académie des Sciences, 73, 147, 237.
[2]	Vrij, A. (1966) Possible Mechanism for the Spontaneous Rupture of Thin, Free Liquid Films. Discussions of the Faraday Society, 42, 23. https://doi.org/10.1039/df9664200023
[3]	Oron, A., Davis, S.H. and Bankoff, S.G. (1997) Long-Scale Evolution of Thin Liquid Films. Reviews of Modern Physics, 69, 931. https://doi.org/10.1103/RevModPhys.69.931
[4]	Bestehorn, M. (2013) Laterally Extended Thin Liquid Films with Inertia under External Vibrations. Physics of Fluids, 25, 114106. https://doi.org/10.1063/1.4830255
[5]	Bestehorn, M. (2009) Fluid Dynamics and Pattern Formation. In: Meyers, R.A., Ed., Contribution to Encyclopedia of Complexity and System Science, 4, Springer, Berlin, 3611. https://doi.org/10.1007/978-0-387-30440-3_214
[6]	Stokes, G.G. (1847) On the Theory of Oscillatory Waves. Trans. Camb. Phil. Soc., 8, 441.
[7]	Dean, R.G. (1965) Stream Function Representation of Non-Linear Ocean Waves. Journal of Geophysical Research, 70, 4561. https://doi.org/10.1029/JZ070i018p04561
[8]	Benjamin, T.B. and Feir, J.E. (1967) The Disintegration of Wave Trains on Deep Water. Journal of Fluid Mechanics, 27, 417. https://doi.org/10.1017/S002211206700045X
[9]	Longuet-Higgins, M.S. and Cokelet, E.D. (1976) The Deformation of Steep Surface Waves on Water. I. A Numerical Method of Computation. Proceedings of the Royal Society of London, A350, 1. https://doi.org/10.1098/rspa.1976.0092
[10]	Zakharov, V.E. (1968) Stability of Periodic Waves of Finite Amplitude on the Surface of a Deep Fluid. Journal of Applied Mechanics and Technical Physics, 9, 190. https://doi.org/10.1007/BF00913182
[11]	Newell, A.C. (1985) Solitons in Mathematics and Physics. CBMS-NSF Regional Conference Series, 48.
[12]	West, B.J., Brueckner, K.A., Janda, R.S., Milder, D.M. and Milton, R.L. (1987) A New Numerical Method for Surface Hydrodynamics. Journal of Geophysical Research, 92, 11803. https://doi.org/10.1029/JC092iC11p11803
[13]	Dommermuth, D.G. and Yue, D. (1987) A High-Order Spectral Method for the Study of Nonlinear Gravity Waves. Journal of Fluid Mechanics, 184, 267. https://doi.org/10.1017/S002211208700288X
[14]	Kharif, C., Pelinovsky, E. and Slunyaev, A. (2009) Rogue Waves in the Ocean. Springer.
[15]	Kwasnicki, M. (2017) Ten Definitions of the Fractional Laplacian Operator. Fractional Calculus and Applied Analysis, 20, 7. https://doi.org/10.1515/fca-2017-0002
[16]	Le Méhauté, B. (1976) An Introduction to Hydrodynamics and Water Waves. Springer, Berlin. https://doi.org/10.1007/978-3-642-85567-2
[17]	Swarztrauber, P.N. (1982) Vectorizing the Fast Fourier Transforms. Academic Press, Cambridge.
[18]	Bestehorn, M. (2018) Computational Physics, with Worked out Examples in FORTRAN and MATLAB, De Gruyter. https://doi.org/10.1515/9783110515145
[19]	Kundu, P.K. and Cohen, I.M. (2004) Fluid Mechanics. Elsevier Academic Press, Cambridge.
[20]	Witham, G.B. (1974) Linear and Nonlinear Waves. John Wiley, Hoboken.
[21]	Zakharov, V.E., Dyachenko, A.I. and Prokofiev, A.O. (2006) Freak Waves as Nonlinear Stage of the Stokes Wave Modulation Instability. European Journal of Mechanics, 25, 677. https://doi.org/10.1016/j.euromechflu.2006.03.004