A simple method of depressing numerical dissipation effects during wave simulation within the Euler model

A simple method of depressing numerical dissipation effects during wave simulation within the Euler model

PDF

Zhe Hu¹, Xiaoying Zhang¹^,^*, Weicheng Cui¹^,²^,³, Fang Wang³, Xiaowen Li¹, Yan Li¹

Acta Oceanologica Sinica | 2020, 39(1) : 141 - 156

Less

Acta Oceanologica Sinica | 2020, 39(1): 141-156

• Marine Technology •

A simple method of depressing numerical dissipation effects during wave simulation within the Euler model

Full

Zhe Hu¹, Xiaoying Zhang¹^,^*, Weicheng Cui¹^,²^,³, Fang Wang³, Xiaowen Li¹, Yan Li¹

Affiliations

¹ Key Laboratory of Ships and Ocean Engineering of Fujian Province, School of Marine Engineering, Jimei University, Xiamen 361021, China

² Deep Sea Technology Research Center, School of Engineering, West Lake University, Hangzhou 310024, China

³ Hadal Science and Technology Research Center, College of Marine Sciences, Shanghai Ocean University, Shanghai 201306, China

Published: 2020-01-25 doi: 10.1007/s13131-019-1524-1

Outline

Abstract

Less

Numerical wave tanks are widely-acknowledged tools in studying waves and wave-structure interactions. They can generate waves under realistic scales and offers more information on the fluid field. However, most numerical wave tanks suffer from issues known as the numerical dissipation and numerical dispersion. The former causes wave energy to be slowly dissipated and the latter shifts wave frequencies during wave propagation. This paper proposes a simple method of depressing numerical dissipation effects on the basis of solving Euler equations using the finite difference method (FDM). The wave propagation solutions are solved analytically taking into account the influence of the damping terms. The main idea of the method is to append a source term to the momentum equation, whose strength is determined by how strong the numerical damping effect is. The method is verified by successfully depressing numerical effects during the simulation of regular linear waves, Stokes waves and irregular waves. By applying the method, wave energy is able to be close to its initial value after long distance of travel.

Key words

numerical dissipation / numerical wave tank / wave simulation / numerical damping reduction / finite difference method

Cite this Article

Zhe Hu, Xiaoying Zhang, Weicheng Cui, Fang Wang, Xiaowen Li, Yan Li. A simple method of depressing numerical dissipation effects during wave simulation within the Euler model[J]. Acta Oceanologica Sinica, 2020 , 39 (1) : 141 -156 . DOI: 10.1007/s13131-019-1524-1

Full Text

Less

1 Introduction

Less

Numerical wave tanks (NWTs) have been widely adopted in studying water waves and wave-related issues (Hasan et al., 2018; Abbasnia et al., 2017; Saincher and Banerjeea, 2015). Compared with analytical methods, NWTs offer approaches to investigating complex phenomena such as high nonlinearity, surface overturning and wave breaking (Anbarsooz et al., 2013; Ma et al., 2016). Being different from experimental methods, NWTs require less costs and offer more detailed information of the fluid field. Within NWTs, tests can be carried out under any conditions, some of which may be too rigorous to be realized in a physical experiment (for instance, Hu et al., 2015).

There are several models which can be frequently found in NWTs for wave simulation. The Euler model and Navier-Stokes (N-S) model may be the most widely used ones, in which the wave motion is described as a combined action of gravity, pressure, inertia and even surface tension, with no viscosity or minor viscosity taken into account (Park et al, 2004; Panicker et al., 2015). When turbulence occurs, which can be found in wave-structure interaction issues, the Reynolds Averaged Navier-Stokes (RANS) model or other models which contain turbulence models may be used (Bihs et al., 2016; Elhanafi et al., 2017). Under the hypothesis of potential flow, the Laplace equation can be adopted and solved with little computational cost (Abbasnia and Soares, 2018; Abbasnia and Ghiasi, 2015). Under certain assumption and special consideration, the general Euler or Laplace model can be simplified into models such as the Boussinesq, Schrödinger and Korteweg-de-Vries (KdV) models. Solving these models requires even lower computational cost, and may gain better results when simulating certain phenomena or revealing some special effects.

Being different from general hydrodynamic problems, some special topics need to be addressed in NWTs, which involve surface capture, wave generation and wave absorption. The volume-of-fluid (VOF) approach is a widely-used surface reconstruction method in NWTs based on solving N-S (or RANS) equations with the finite difference method (FDM) or finite volume method (FVM). Other surface processing techniques contain the level set method (Schillaci et al., 2016), marker and cell (MAC) method (De Paulo et al., 2007) and particle method such as smooth particle hydrodynamics (SPH) method (Alvarado-Rodríguez et al., 2017), etc. Wave generation is usually realized by setting the boundary condition of the flow field in a pre-designed manner, or by appending source terms to conservation equations within certain regions. Typically, the imitation of physical wave-generating paddles and the free income boundary (Hu et al., 2017) are wave makers of the first category, yet the mass and momentum source method are ones of the second category (Liu et al., 2015). Similarly, wave absorption approaches can be also divided into two categories, namely by setting boundary conditions (e.g., the open boundary method) or making use of source terms (e.g., the imitation of spongy layer or porous media). It should be mentioned that the wave absorption method had better match wave maker in one NWT and vice versa. For instance, if the free income boundary is adopted for wave generation, then absorption methods which can discharge fluid mass (e.g., the open boundary method) are preferred, as the free income boundary will induce a slow transportation of fluid mass and therefore change the mean free-surface level after long-term simulation (Dean and Dalrymple, 1991).

Based on above statements, a NWT can be established only after the above-mentioned techniques are well prepared, together with a reliable flow solver. However, even if one confirms that all parts of the NWT act correctly, one may still find it difficult to simulate certain waves within the NWT. One important reason leading to this difficulty is the so-called numerical dissipation and dispersion. It is caused by the fact that the numerical solution usually exhibits deviation from the real solution of the partial differences equations (PDEs), and most of the time the resulting error is dominated by either dissipative or dispersive features. During wave simulation, numerical dissipation makes the numerical wave amplitude smaller than the intended target one, and numerical dispersion changes the dispersion relationship of the wave. Both effects should be avoided in the study of wave propagation mechanisms, as they may lead to non-physical phenomena or give misleading information on wave features. Whereas the tricky thing is that numerical dissipation and dispersion can only be declined or depressed, but can be never completely eliminated.

To eliminate the effect of numerical dissipation, various techniques have been proposed in literature. Nazari et al. (2015) proposed a high-order diagonally implicit Runge-Kutta scheme with low dispersion and negligible dissipation. Li et al. (2015) proposed a low dissipation numerical scheme for Implicit Large Eddy Simulation by combining two schemes for reconstruction and the Riemann solver respectively. Schranner et al. (2015) proposed a method for quantifying the effective numerical dissipation rate and effective numerical viscosity in Computational Fluid Dynamics (CFD) simulations, which allowed the determination of numerical dissipation rates locally in a physical-space representation. Cao et al. (2017) invented a dissipation-free numerical method to capture shock wave propagation in 1-dimensional fluid flow problems, which achieved a higher accuracy while using much less grid number. Beljadid et al. (2017) proposed a class of numerical schemes for the approximation of weak solutions to nonlinear hyperbolic systems in nonconservative form with well-controlled dissipation. Soares (2019) proposed a new explicit-implicit time integration technique for wave propagation analysis in which controllable algorithm dissipation is provided.

This paper investigates the effect of numerical dissipation upon wave simulation in NWTs under the framework of the Euler equations and further gives a simple method of depressing the influence of numerical dissipation during wave propagation. Conventional numerical dissipation combating methods, as has been mentioned above, tend to use special schemes or new numerical algorithms. By contrast, in our method a source term is appended to the momentum equation which compensates for the energy loss due to numerical dissipation, and hence there is no need to modify the original numerical schemes. Another difference is that our method is specially designed for simulating water surface waves, which can be easily implemented and requires low computational cost.

In this paper, the corrected PDEs of the adopted numerical algorithm are derived and simplified. Wave-type solutions of the corrected PDEs are deduced and analyzed. On the basis of the wave solution, a simple method is proposed to reduce the energy dissipation induced by numerical damping by means of appending enhancing source terms to the momentum equations. The method is verified by simulating regular linear waves, nonlinear Stokes waves and irregular wave spectrums. The energy of the simulated waves barely changes after long distance of propagation, if the method is adopted.

2 Corrected PDEs of the numerical scheme

Less

Most of the time, the numerical solutions of PDEs present certain deviation from the PDEs’ precise solutions. In other words, the equation system represented by the numerical method, also called the corrected PDEs, is somewhat different from the original PDEs. Unless for special purposes, the corrected PDEs should be close to the original PDEs when the size of the division (or discretization) of the numerical method becomes small. Usually, determining the corrected PDEs of the numerical method is the first step of evaluating numerical dissipation. Since different numerical methods correspond to different corrected PDEs, without loss of generality, a simple FDM algorithm with staggered grids and even mesh is adopted for analysis.

Standard Euler equations within two-dimensional space can be written as follows.

(1a)

$ \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} = 0,$

(1b)

$\tag{1b}\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial z}} = - \frac{1}{\rho }\frac{{\partial P}}{{\partial x}},$

(1c)

$\tag{1c}\frac{{\partial v}}{{\partial t}} + u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial z}} = - \frac{1}{\rho }\frac{{\partial P}}{{\partial z}} - g,$

where

$u$

and

$v$

are the x- and z-directional component of the velocity vector,

$P$

is pressure,

$g$

is the gravitational acceleration, and

$\rho $

is the fluid density. As is shown in Fig. 1, staggered grid is used where velocities are defined at the centers of cell boundaries and pressures at cell centers. Under the notation of Fig. 1, Eq. (1b) and Eq. (1c) are spatially discretized at position (i+1/2, j) and (i, j+1/2), respectively. Yet the numerical discretization of Eq. (1a) is performed at position (i, j).

In this way, discretizing the convective terms of Eq. (1b) and Eq. (1c) with central difference scheme (CDS) yields

(2a)

$\tag{2a}{\bf{D}}{\left({u\frac{{\partial u}}{{\partial x}}} \right)_{i + 1/2,j}} = {u_{i + 1/2,j}}\frac{{{u_{i + 3/2,j}} - {u_{i - 1/2,j}}}}{{2\Delta x}},$

(2b)

$\tag{2b} \begin{split} {\bf{D}}{\left( {v\frac{{\partial u}}{{\partial z}}} \right)_{i + 1/2,j}} = & \frac{{{v_{i + 1,j + 1/2}} + {v_{i + 1,j - 1/2}} + {v_{i,j + 1/2}} + {v_{i,j - 1/2}}}}{4} \times\\ & \frac{{{u_{i + 1/2,j + 1}} - {u_{i + 1/2,j - 1}}}}{{2\Delta z}}, \end{split} $

(2c)

$\tag{2c} \begin{split} {\bf{D}}{\left( {u\frac{{\partial v}}{{\partial x}}} \right)_{i,j + 1/2}} =& \frac{{{u_{i + 1/2,j + 1}} + {u_{i - 1/2,j + 1}} + {u_{i + 1/2,j}} + {u_{i - 1/2,j}}}}{4} \times\\ & \frac{{{v_{i + 1,j + 1/2}} - {v_{i - 1,j + 1/2}}}}{{2\Delta x}}, \end{split} $

(2d)

$\tag{2d}{\bf{D}}{\left({v\frac{{\partial v}}{{\partial z}}} \right)_{i,j + 1/2}} = {v_{i,j + 1/2}}\frac{{{v_{i,j + 3/2}} - {v_{i,j - 1/2}}}}{{2\Delta z}},$

where

$\Delta x$

and

$\Delta z$

are the x- and z-directional grid size, and operator

${\bf{D}}\left({\;} \right)$

denotes performing (spatial or temporal) numerical discretization upon the bracketed terms. By performing Taylor expansions (Eq. (3) for example) with respect to the intended locations, Eq. (2) can be converted into the summations of certain series, which are formulated in Eq. (4).

(3)

$ \begin{split} {u_{i - 1/2,j}} =& {u_{i + 1/2,j}} - {\left( {\frac{{\partial u}}{{\partial x}}} \right)_{i + 1/2,j}}\Delta x + \frac{1}{2}{\left( {\frac{{\partial u}}{{\partial x}}} \right)^2}_{i + 1/2,j}\Delta {x^2} + ... +\\ & \frac{{{{\left({ - 1} \right)}^n}}}{{n!}}{\left({\frac{{\partial u}}{{\partial x}}} \right)^n}_{i + 1/2,j}\;\;\Delta {x^n} + O\left({\Delta {x^{n + 1}}} \right). \end{split} $

(4a)

$\tag{4a}{\bf{D}}{\left({u\frac{{\partial u}}{{\partial x}}} \right)_{i + 1/2,j}} = {\left. {u\left({\frac{{\partial u}}{{\partial x}} + \sum\limits_{n = 3,5,7, \cdots } {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {x^n}}}\Delta {x^{n - 1}}} } \right)} \right|_{i + 1/2,j}},$

(4b)

$\tag{4b}\begin{split} & {\bf{D}}{\left({v\frac{{\partial u}}{{\partial z}}} \right)_{i + 1/2,j}} = \\& \left({v + \sum\limits_{n = 2,4,6, \cdots } {\frac{1}{{{2^n}n!}}\sum\limits_{k = 0,2,4, \cdots }^n {C_n^k\Delta {x^k}\Delta {z^{n - k}}} \frac{{{\partial ^n}v}}{{\partial {x^k}\partial {z^{n - k}}}}} } \right) \times\\ & \left.\left({\frac{{\partial u}}{{\partial z}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {z^n}}}\Delta {z^{n - 1}}} } \right) \right|_{i + 1/2,j},\end{split}$

(4c)

$\tag{4c}\begin{split} & {\bf{D}}{\left({u\frac{{\partial v}}{{\partial y}}} \right)_{i,j + 1/2}} = \\ & \left({u + \sum\limits_{n = 2,4,6, \cdots } {\frac{1}{{{2^n}n!}}\sum\limits_{k = 0,2,4, \cdots }^n {C_n^k\Delta {x^k}\Delta {z^{n - k}}} \frac{{{\partial ^n}u}}{{\partial {x^k}\partial {z^{n - k}}}}} } \right) \times\\& \left. \left({\frac{{\partial v}}{{\partial x}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {x^n}}}\Delta {x^{n - 1}}} } \right) \right|_{i,j + 1/2} ,\end{split}$

(4d)

$\tag{4d}{\bf{D}}{\left({v\frac{{\partial v}}{{\partial z}}} \right)_{i,j + 1/2}} = {\left. {v\left({\frac{{\partial v}}{{\partial z}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {z^n}}}\Delta {z^{n - 1}}} } \right)} \right|_{i,j + 1/2}}.$

The discretization of the pressure gradient and velocity divergence also follows the CDS, which can be written as

(5a)

$\tag{5a} \begin{split} {\bf{D}}{\left({\frac{{\partial P}}{{\partial x}}} \right)_{i + 1/2,j}} =& \frac{{{P_{i + 1,j}} - {P_{i,j}}}}{{\Delta x}} \\ =& {\left. {\left[ {\frac{{\partial P}}{{\partial x}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}P}}{{\partial {x^n}}}{{\left({\frac{{\Delta x}}{2}} \right)}^{n - 1}}} } \right]} \right|_{i + 1/2,j}}, \end{split} $

(5b)

$\tag{5b} \begin{split} {\bf{D}}{\left({\frac{{\partial P}}{{\partial z}}} \right)_{i,j + 1/2}} =& \frac{{{P_{i,j + 1}} - {P_{i,j}}}}{{\Delta z}} \\ =& {\left. {\left[ {\frac{{\partial P}}{{\partial z}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}P}}{{\partial {z^n}}}{{\left({\frac{{\Delta z}}{2}} \right)}^{n - 1}}} } \right]} \right|_{i,j + 1/2}}, \end{split} $

(5c)

$\tag{5c}\begin{split} {\bf{D}}{\left({\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}}} \right)_{i,j}} =& \frac{{{u_{i + 1/2,j}} - {u_{i - 1/2,j}}}}{{\Delta x}} + \frac{{{v_{i,j + 1/2}} - {v_{i,j - 1/2}}}}{{\Delta z}} \\ =& \left[ \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} + \sum\limits_{n = 3,5,7, \cdots } {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {x^n}}}{{\left({\frac{{\Delta x}}{2}} \right)}^{n - 1}}} +\right.\\& \left.\left.\sum\limits_{n = 3,5,7, \cdots } {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {z^n}}}{{\left({\frac{{\Delta z}}{2}} \right)}^{n - 1}}} \right] \right|_{i,j}.\end{split} $

A simple forward Euler scheme is used for time advancing, and can be written as

(6)

${\bf{D}}{\left({\frac{{\partial u}}{{\partial t}}} \right)^t} = \frac{{{u^{t + \Delta t}} - {u^t}}}{{\Delta t}} = {\left. {\left({\frac{{\partial u}}{{\partial t}} + \sum\limits_{n = 2}^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {t^n}}}\Delta {t^{n - 1}}} } \right)} \right|^t}.$

The coupled velocity-pressure equations can be solved using a projection method, in which the momentum equations can be formulated as

(7a)

$\tag{7a}\begin{split} {\bf{D}}\left({\frac{{\partial u}}{{\partial t}}} \right)_{i + 1/2,j}^t =& - \frac{1}{\rho }{\bf{D}}\left({\frac{{\partial P}}{{\partial x}}} \right)_{i + 1/2,j}^{t + \Delta t} -\\ & {\bf{D}}{\left({u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial z}}} \right)^t}_{i + 1/2,j}, \end{split} $

(7b)

$\tag{7b} \begin{split} {\bf{D}}\left({\frac{{\partial v}}{{\partial t}}} \right)_{i,j + 1/2}^t =& - \frac{1}{\rho }{\bf{D}}\left({\frac{{\partial P}}{{\partial z}}} \right)_{i,j + 1/2}^{t + \Delta t} -\\ & {\bf{D}}{\left({u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial z}}} \right)^t}_{i,j + 1/2} - g. \end{split} $

Noticing that the pressure gradient term is performed at time step

$t + \Delta t$

, another Taylor expansion with respect to t is required to keep all terms temporally consistent, and the resulted corrected momentum equations are

(8a)

$\tag{8a} \begin{split} \frac{{\partial u}}{{\partial t}} +& u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial z}} + \frac{{\partial P}}{{\partial x}} = - \sum\limits_{n = 2}^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {t^n}}}\Delta {t^{n - 1}}} -\\& \left\{ \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}P}}{{\partial {x^n}}}{{\left({\frac{{\Delta x}}{2}} \right)}^{n - 1}}} + \sum\limits_{m = 1}^\infty \frac{1}{{m!}}\frac{{{\partial ^m}}}{{\partial {t^m}}}\right.\times\\& \left.\left[ {\frac{{\partial P}}{{\partial x}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}P}}{{\partial {x^n}}}{{\left({\frac{{\Delta x}}{2}} \right)}^{n - 1}}} } \right]\Delta {t^m} \right\}- \\& u\sum\limits_{n = 3,5,7, \cdots } {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {x^n}}}\Delta {x^{n - 1}}} - v\sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {z^n}}}\Delta {z^{n - 1}}}- \\ & \frac{{\partial u}}{{\partial z}}\sum\limits_{n = 2,4,6, \cdots } {\frac{1}{{{2^n}n!}}\sum\limits_{k = 0,2,4, \cdots }^n {C_n^k\Delta {x^k}\Delta {z^{n - k}}} \frac{{{\partial ^n}v}}{{\partial {x^k}\partial {z^{n - k}}}}} -\\ & \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}u}}{{\partial {z^n}}}\Delta {z^{n - 1}}} \sum\limits_{n = 2,4,6, \cdots } \frac{1}{{{2^n}n!}}\times \\ &{\sum\limits_{k = 0,2,4, \cdots }^n {C_n^k\Delta {x^k}\Delta {z^{n - k}}} \frac{{{\partial ^n}v}}{{\partial {x^k}\partial {z^{n - k}}}}},\end{split}$

(8b)

$\tag{8b}\begin{split} \frac{{\partial v}}{{\partial t}} +& u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial z}} + \frac{{\partial P}}{{\partial z}} + g = - \sum\limits_{n = 2}^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {t^n}}}\Delta {t^{n - 1}}} -\\ & \left\{ \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}P}}{{\partial {z^n}}}{{\left({\frac{{\Delta z}}{2}} \right)}^{n - 1}}} + \sum\limits_{m = 1}^\infty \frac{1}{{m!}}\right.\times\\ & \left.\frac{{{\partial ^m}}}{{\partial {t^m}}}\left[ {\frac{{\partial P}}{{\partial z}} + \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}P}}{{\partial {z^n}}}{{\left({\frac{{\Delta z}}{2}} \right)}^{n - 1}}} } \right]\Delta {t^m} \right\}- \\& v\sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {z^n}}}\Delta {z^{n - 1}}} - u\sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {x^n}}}\Delta {x^{n - 1}}} -\\& \frac{{\partial v}}{{\partial x}}\sum\limits_{n = 2,4,6, \cdots } {\frac{1}{{{2^n}n!}}\sum\limits_{k = 0,2,4, \cdots }^n {C_n^k\Delta {x^k}\Delta {z^{n - k}}} \frac{{{\partial ^n}u}}{{\partial {x^k}\partial {z^{n - k}}}}} -\\& \sum\limits_{n = 3,5,7, \cdots }^\infty {\frac{1}{{n!}}\frac{{{\partial ^n}v}}{{\partial {x^n}}}\Delta {x^{n - 1}}} \sum\limits_{n = 2,4,6, \cdots }\frac{1}{{{2^n}n!}}\times\\& {\sum\limits_{k = 0,2,4, \cdots }^n {C_n^k\Delta {x^k}\Delta {z^{n - k}}} \frac{{{\partial ^n}u}}{{\partial {x^k}\partial {z^{n - k}}}}} .\end{split} $

The corrected continuity equation can be formulated as

(8c)

$\tag{8c} \begin{split} \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} =& - \sum\limits_{n = 3,5,7, \cdots } {\frac{1}{{n!}}\left({\frac{{{\partial ^n}u}}{{\partial {x^n}}}} \right){{\left({\frac{{\Delta x}}{2}} \right)}^{n - 1}}} -\\ & \sum\limits_{n = 3,5,7, \cdots } {\frac{1}{{n!}}\left({\frac{{{\partial ^n}v}}{{\partial {z^n}}}} \right){{\left({\frac{{\Delta z}}{2}} \right)}^{n - 1}}} . \end{split} $

By preserving only the leading-order terms of the right side of Eq. (8), the corrected PDEs can be reformulated as

(9a)

$\tag{9a} \begin{split} \frac{{\partial u}}{{\partial t}} +& u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial z}} + \frac{{\partial P}}{{\partial x}} = - \frac{{\Delta t}}{2}\frac{{{\partial ^2}u}}{{\partial {t^2}}} - \frac{{\Delta {x^2}}}{{24}}\frac{{{\partial ^3}P}}{{\partial {x^3}}} -\\& \Delta t\frac{{{\partial ^2}P}}{{\partial x\partial t}} - u\frac{{\Delta {x^2}}}{6}\frac{{{\partial ^3}u}}{{\partial {x^3}}} - v\frac{{\Delta {z^2}}}{6}\frac{{{\partial ^3}u}}{{\partial {z^3}}}- \\& \frac{1}{8}\frac{{\partial u}}{{\partial z}}\left({\Delta {z^2}\frac{{{\partial ^2}v}}{{\partial {z^2}}} + \Delta {x^2}\frac{{{\partial ^2}v}}{{\partial {x^2}}}} \right) +\\& O\left({\Delta {t^2},\Delta t\Delta {x^2},\Delta {x^4},\Delta {z^4},\Delta {x^2}\Delta {z^2}} \right),\end{split} $

(9b)

$ \tag{9b} \begin{split} \frac{{\partial v}}{{\partial t}} +& u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial z}} + \frac{{\partial P}}{{\partial z}} + g = - \frac{{\Delta t}}{2}\frac{{{\partial ^2}v}}{{\partial {t^2}}} - \frac{{\Delta {z^2}}}{{24}}\frac{{{\partial ^3}P}}{{\partial {z^3}}} -\\& \Delta t\frac{{{\partial ^2}P}}{{\partial z\partial t}} - v\frac{{\Delta {z^2}}}{6}\frac{{{\partial ^3}v}}{{\partial {z^3}}} - u\frac{{\Delta {x^2}}}{6}\frac{{{\partial ^3}v}}{{\partial {x^3}}} -\\& \frac{1}{8}\frac{{\partial v}}{{\partial x}}\left({\Delta {z^2}\frac{{{\partial ^2}u}}{{\partial {z^2}}} + \Delta {x^2}\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right) +\\& O\left({\Delta {t^2},\Delta t\Delta {z^2},\Delta {x^4},\Delta {z^4},\Delta {x^2}\Delta {z^2}} \right) ,\end{split} $

(9c)

$\tag{9c}\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} = - \frac{{\Delta {x^2}}}{{24}}\frac{{{\partial ^3}u}}{{\partial {x^3}}} - \frac{{\Delta {z^2}}}{{24}}\frac{{{\partial ^3}v}}{{\partial {z^3}}} + O\left({\Delta {x^4},\Delta {z^4}} \right).$

3 Wave propagation and a simple method of depressing dissipation

Less

In numerical methods, the systems to be solved are the corrected PDEs rather than the original PDEs. Like traditional PDEs, the solutions of the corrected PDEs are not unique, and the influence of numerical dissipation may be multiple. This session discusses how wave solution is changed with the existence of numerical dissipation.

3.1 Further simplification upon the corrected PDEs

Let us consider a small-amplitude regular wave, which is derived on the basis of the hypothesis of small velocity and dynamic pressure. Under such assumption, nonlinear terms of Euler equations are considered to be of higher orders, and hence Eq. (9) can be simplified as

(10a)

$\tag{10a} \begin{split} \frac{{\partial u}}{{\partial t}} + \frac{{\partial P}}{{\partial x}} =& - \frac{{\Delta t}}{2}\frac{{{\partial ^2}u}}{{\partial {t^2}}} - \frac{{\Delta {x^2}}}{{24}}\frac{{{\partial ^3}P}}{{\partial {x^3}}} - \Delta t\frac{{{\partial ^2}P}}{{\partial x\partial t}} + \\ & O\left({\Delta {t^2},\Delta t\Delta {x^2},\Delta {x^4},\Delta {z^4},\Delta {x^2}\Delta {z^2}} \right), \end{split} $

(10b)

$\tag{10b} \begin{split} \frac{{\partial v}}{{\partial t}} + \frac{{\partial P}}{{\partial z}} + g =& - \frac{{\Delta t}}{2}\frac{{{\partial ^2}v}}{{\partial {t^2}}} - \frac{{\Delta {z^2}}}{{24}}\frac{{{\partial ^3}P}}{{\partial {z^3}}} - \Delta t\frac{{{\partial ^2}P}}{{\partial z\partial t}} +\\ & O\left({\Delta {t^2},\Delta t\Delta {z^2},\Delta {x^4},\Delta {z^4},\Delta {x^2}\Delta {z^2}} \right), \end{split} $

(10c)

$\tag{10c}\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} = - \frac{{\Delta {x^2}}}{{24}}\frac{{{\partial ^3}u}}{{\partial {x^3}}} - \frac{{\Delta {z^2}}}{{24}}\frac{{{\partial ^3}v}}{{\partial {z^3}}} + O\left({\Delta {x^4},\Delta {z^4}} \right).$

During wave simulation, the grid size and time interval should be small enough compared with the wave length and wave period. This fact can be represented as

$\omega \Delta t \sim O\left(\varepsilon \right)$

and

$k\Delta x\;(.or.\;k\Delta z) \sim O\left(\varepsilon \right)$

, where

$\varepsilon $

is a small parameter. Using this fact, the terms on the right sides of Eq. (10) have following orders

(11a)

$\tag{11a}\frac{{\Delta t}}{2}\frac{{{\partial ^2}u}}{{\partial {t^2}}}\; \sim \varepsilon \frac{{\partial u}}{{\partial t}},\;\;\;\;\frac{{\rm \Delta t}}{2}\frac{{{\partial ^2}v}}{{\partial {t^2}}}\; \sim \varepsilon \frac{{\partial v}}{{\partial t}},$

(11b)

$\tag{11b}\frac{{\Delta {x^2}}}{{24}}\frac{{{\partial ^3}P}}{{\partial {x^3}}} \sim \;{\varepsilon ^2}\frac{{\partial P}}{{\partial x}},\;\;\;\;\frac{{\Delta {z^2}}}{{24}}\frac{{{\partial ^3}P}}{{\partial {z^3}}}\; \sim {\varepsilon ^2}\frac{{\partial P}}{{\partial z}},$

(11c)

$\tag{11c}\Delta t\frac{{{\partial ^2}P}}{{\partial x\partial t}}\; \sim \;\varepsilon \frac{{\partial P}}{{\partial x}},\;\;\;\;\rm\Delta t\frac{{{\partial ^2}P}}{{\partial z\partial t}} \sim \varepsilon \frac{{\partial P}}{{\partial z}},$

(11d)

$\tag{11d}\frac{{\Delta {x^2}}}{{24}}\frac{{{\partial ^3}u}}{{\partial {x^3}}} \sim \;{\varepsilon ^2}\frac{{\partial u}}{{\partial x}},\;\;\;\;\;\;\;\frac{{\Delta {z^2}}}{{24}}\frac{{{\partial ^3}v}}{{\partial {z^3}}} \sim {\varepsilon ^2}\frac{{\partial v}}{{\partial z}}.$

Omitting terms of order

${\varepsilon ^2}$

, Eq. (10) can be converted into the following equations.

(12a)

$\tag{12a}\frac{{\partial u}}{{\partial t}} + \frac{{\partial P}}{{\partial x}} = - \frac{{\Delta t}}{2}\frac{{{\partial ^2}u}}{{\partial {t^2}}} - \Delta t\frac{{{\partial ^2}P}}{{\partial x\partial t}},$

(12b)

$\tag{12b}\frac{{\partial v}}{{\partial t}} + \frac{{\partial P}}{{\partial z}} + g = - \frac{{\Delta t}}{2}\frac{{{\partial ^2}v}}{{\partial {t^2}}} - \Delta t\frac{{{\partial ^2}P}}{{\partial z\partial t}},$

(12c)

$\tag{12c}\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} = 0.$

The second terms of the right sides of Eq. (12a) and Eq. (12b) can be further transformed as follows:

(13a)

$\tag{13a} \begin{split} - \Delta t\frac{{{\partial ^2}P}}{{\partial x\partial t}} = & - \Delta t\frac{\partial }{{\partial t}}\left({ - \frac{{\partial u}}{{\partial t}} - \frac{{\Delta t}}{2}\frac{{{\partial ^2}u}}{{\partial {t^2}}} - \Delta t\frac{{{\partial ^2}P}}{{\partial x\partial t}}} \right) \\ =& \Delta t\frac{{{\partial ^2}u}}{{\partial {t^2}}} + O\left({{\varepsilon ^2}} \right), \end{split} $

(13a)

$\tag{13a} \begin{split} - \Delta t\frac{{{\partial ^2}P}}{{\partial z\partial t}} =& - \Delta t\frac{\partial }{{\partial t}}\left({ - \frac{{\partial v}}{{\partial t}} - \frac{{\Delta t}}{2}\frac{{{\partial ^2}v}}{{\partial {t^2}}} - \Delta t\frac{{{\partial ^2}P}}{{\partial z\partial t}}} \right) \\ =& \Delta t\frac{{{\partial ^2}v}}{{\partial {t^2}}} + O\left({{\varepsilon ^2}} \right). \end{split} $

Noticing the right sides of Eq. (12) only preserve terms of order

$\varepsilon $

, the corrected PDEs are finally rewritten as

(14a)

$\tag{14a}\frac{{\partial u}}{{\partial t}} + \frac{{\partial P}}{{\partial x}} = \frac{{\Delta t}}{2}\frac{{{\partial ^2}u}}{{\partial {t^2}}},$

(14b)

$\tag{14b}\frac{{\partial v}}{{\partial t}} + \frac{{\partial P}}{{\partial z}} + g = \frac{{\Delta t}}{2}\frac{{{\partial ^2}v}}{{\partial {t^2}}},$

(14c)

$\tag{14c}\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial z}} = 0.$

The corrected PDEs differ from the original PDEs in terms on the right side of Eq. (14), which will be proved to be numerical dissipative later.

3.2 Wave solutions

Under the hypothesis of being irrotational, Eq. (14) can be expressed with the velocity potential as

(15)

$\left\{ \begin{aligned} & {\nabla ^2}\Phi = 0 \\ & \frac{{\partial \Phi }}{{\partial t}} + gz + P - \frac{{\Delta t}}{2}\frac{{{\partial ^2}\Phi }}{{\partial {t^2}}} = 0\end{aligned} \right.,$

where

$\Phi $

is the velocity potential. Taking into account the kinematic free surface condition and bottom condition, the governing equations can be written as

(16)

$\left\{ \begin{aligned} & {\nabla ^2}\Phi = 0\quad\quad\quad\quad - h < z < 0 \\& \frac{{\partial \Phi }}{{\partial t}} + g\eta - \frac{{\Delta t}}{2}\frac{{{\partial ^2}\Phi }}{{\partial {t^2}}} = 0\quad\quad\quad\quad z = 0 \\& \frac{{\partial \eta }}{{\partial t}} = \frac{{\partial \Phi }}{{\partial z}}\quad\quad\quad\quad\quad\quad z = 0 \\& \frac{{\partial \Phi }}{{\partial z}} = 0\quad\quad\quad\quad z = - h \end{aligned} \right.,$

where h is water depth. The unified free surface condition can be obtained by combining the kinematic free surface condition and the dynamic free surface condition, and is written as

(17)

$\frac{{{\partial ^2}\Phi }}{{\partial {t^2}}} + g\frac{{\partial \Phi }}{{\partial z}} - \frac{{\Delta t}}{2}\frac{{{\partial ^3}\Phi }}{{\partial {t^3}}} = 0\qquad\quad z = 0.$

The wave-type solution of Eq. (16) in complex form is

(18a)

$\tag{18a}\Phi = - a\frac{{{\rm i}g}}{\omega }\frac{{\cosh k\left({z + h} \right)}}{{\cosh kh}}{{\rm e}^{{\rm i}\left({kx - \omega t} \right)}},$

(18b)

$\tag{18b}\eta = a{{\rm e}^{{\rm i}\left({kx - \omega t} \right)}},$

where i is the complex unit, a is the wave amplitude, k and

$\omega $

are wave number and angular frequency, respectively. Substituting Eq. (18) into Eq. (17), the dispersion equation can be derived as

(19)

$ - \frac{{{\rm i}\Delta t}}{2}{\omega ^3} - {\omega ^2} + gk\tanh kh = 0.$

Equation (19) differs from the conventional dispersion relationship in having the term

$ - {\rm i}\Delta t{\omega ^3}/2$

. Two categories of solutions which contain physical meanings can be drawn from Eq. (19).

(1) Solution with real wave number

$k_0$

If the wave number is real (i.e.,

$k = {k_0}$

), with the help of symbolic operation tools,

$\omega $

can be gained by solving a cubic algebraic equation as

(20a)

$\tag{20a}{\omega _1} = - \frac{1}{{3\Delta t}}\left[ {\sqrt 3 \left({\alpha - \frac{1}{\alpha }} \right) + \left({\alpha + \frac{1}{\alpha } - 2} \right){\rm i}} \right],$

(20b)

$\tag{20b}{\omega _2} = - \frac{1}{{3\Delta t}}\left[ { - \sqrt 3 \left({\alpha - \frac{1}{\alpha }} \right) + \left({\alpha + \frac{1}{\alpha } - 2} \right){\rm i}} \right],$

(20c)

$\tag{20c}{\omega _3} = \frac{2}{{3\Delta t}}\left({\alpha + \frac{1}{\alpha } + 1} \right){\rm i},$

where

(21a)

$\tag{21a}\alpha = {\left[ {1 + \frac{1}{8}\left({27{\omega _0}^2\Delta {t^2} + 3{\omega _0}\Delta t\sqrt {81{\omega _0}^2\Delta {t^2} + 48} } \right)} \right]^{1/3}},$

(21b)

$\tag{21b}{\omega _0}^2 = g{k_0}\tanh {k_0}h.$

${\omega _1}$

and

${\omega _2}$

represent left-traveling and right-traveling waves, respectively.

${\omega _3}$

represents energy growing unboundedly within the field, and thereby should be abandoned.

One may find that Eq. (20) and Eq. (21) are lengthy, and therefore are not convenient for analysis. Noticing the fact that

$ - {\rm i}\Delta t{\omega ^3}/2$

is a small number,

$\omega $

must be closed to

${\omega _0}$

. Let

$\omega = {\omega _0} + \delta $

, and Eq. (19) is rewritten as

(22)

$ - \frac{{{\rm i}\Delta t}}{2}{\left({{\omega _0} + \delta } \right)^3} - {\left({{\omega _0} + \delta } \right)^2} + g{k_0}\tanh {k_0}h = 0.$

Keeping terms corresponding to order

$O\left(\delta \right)$

, the solution of Eq. (22) is

(23)

$\delta = - \frac{{{\rm i}\Delta t{\omega _0}^2}}{{4 + 3{\rm i}\Delta t{\omega _0}}} \to - \frac{{{\rm i}\Delta t{\omega _0}^2}}{4}.$

By applying Eq. (23), the wave solution of the corrected PDEs is written as

(24a)

$\tag{24a}\Phi = - a\frac{{{\rm i}g}}{{{\omega _0}}}{{\rm e}^{ - \;\frac{{\Delta t{\omega _0}^2}}{4}t}}\frac{{\cosh {k_0}\left({z + h} \right)}}{{\cosh {k_0}h}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}},$

(24b)

$\tag{24b}\eta = a{{\rm e}^{ - \;\frac{{\rm\Delta t{\omega _0}^2}}{4}t}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}},$

where

${\omega _0}$

and

${k_0}$

satisfy the conventional dispersion relationship. Eq. (24) indicates that a free-traveling linear monochromatic wave will dissipate exponentially with respect to time, if no extra energy is supplied. The dissipation rate is related to the selected time step and wave frequency. In NWTs, wave makers (e.g., a wave generating paddle) usually work in predesigned manners which are not affected by the inner flow field. Namely, wave makers of NWTs act as energy suppliers to the computational domain. This phenomenon can be explained by solution of the second category.

(2) Solution with real angular frequency

${\omega _0}$

If the angular frequency is real, namely

$\omega = {\omega _0}$

, then the wave number must have a small deviation from

${k_0}$

. Let

$k = {k_0} + \delta $

, and the dispersion equation is written as

(25)

$ - \frac{{{\rm i}\Delta t}}{2}{\omega _0}^3 - {\omega _0}^2 + g\left({{k_0} + \delta } \right)\tanh \left({{k_0} + \delta } \right)h = 0.$

Using Taylor expansion, the third term on the left side of Eq. (25) can be written as

(26)

$\begin{split} g\left({{k_0} + \delta } \right)\tanh \left({{k_0} + \delta } \right)h =& g{k_0}\tanh {k_0}h + \\ &\frac{\partial }{{\partial {k_0}}}\left({g{k_0}\tanh {k_0}h} \right) \cdot \delta + O\left({{\delta ^2}} \right) \\ =& \omega _0^2 + 2\omega _0^{}{C_g} \cdot \delta + O\left({{\delta ^2}} \right), \end{split}$

where

${C_{\rm g}}$

is the group velocity. Substituting Eq. (26) into Eq. (25) yields

(27)

$\delta = \frac{{{\rm i}\Delta t{\omega _0}^2}}{{4{C_g}}}.$

The expressions of velocity potential and surface elevation are finally written as

(28a)

$\tag{28a}\Phi = - a\frac{{{\rm i}g}}{{{\omega _0}}}{{\rm e}^{ - \;\frac{{\Delta t{\omega _0}^2}}{{4{C_{\rm g}}}}x}}\frac{{\cosh {k_0}\left({z + h} \right)}}{{\cosh {k_0}h}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}},$

(28b)

$\tag{28b}\eta = a{{\rm e}^{ - \;\frac{{\rm\Delta t{\omega _0}^2}}{{4{C_{\rm g}}}}x}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}}.$

Equation (28) reveals that the amplitude (energy) of a free-traveling linear monochromatic wave declines exponentially along its propagating direction. In NWTs, wave makers keep donating energy to their neighboring fluid particles, and therefore make waves adjacent to wave makers present relatively large amplitudes. However, as waves are travelling away from wave generators, their energy is reduced by numerical dissipation effects in a manner of Eq. (28). It is worthy of noting that the term

$\dfrac{{\cosh {k_0}\left({z + h} \right)}}{{\cosh {k_0}h}}$

in Eq. (28a) should technically be

$\dfrac{{\cosh \left[ {\left({{k_0} + \delta } \right)\left({z + h} \right)} \right]}}{{\cosh \left[ {\left({{k_0} + \delta } \right)h} \right]}}$

. Yet it can be easily proofed that

$\dfrac{{\cosh \left[ {\left({{k_0} + \delta } \right)\left({z + h} \right)} \right]}}{{\cosh \left[ {\left({{k_0} + \delta } \right)h} \right]}} \to$

$ \dfrac{{\cosh {k_0}\left({z + h} \right)}}{{\cosh {k_0}h}} + O\left({{\delta ^2}} \right) $

under both deep and shallow water depth, and therefore only the leading term is reserved in Eq. (28a).

3.3 A simple method of depressing energy dissipation

The dominate reason which leads to the troubling numerical damping effects is the additional terms on the right side of Eq. (14). For the purpose of declining numerical dissipation (namely, declining the right side of Eq. (14)), the simplest way is to use a denser mesh and a shorter time interval. In addition, one can also try to use high-order differencing schemes. However, both means require high computational costs and more computer memory. Besides, high-order differencing schemes can increase the risk of numerical instability.

Instead of using high-order differencing schemes or small spatial and temporal division upon the solution region, an alternative method is adding additional terms to the original Euler equations (i.e., Eq. (1)), which counteract the damping effect induced by numerical dissipation. These additional terms had better contain no derivatives, since realizing derivatives numerically will again cause additional computational cost and numerical dissipation or dispersion, which is much in common with the method of using high-order schemes.

Let us consider a simple source term

$c{\boldsymbol{u}}$

, where

$c$

is a small constant and

${\boldsymbol{u}}$

is the velocity vector. By involving this source term, the Euler equations are modified as

(29)

$\frac{{\partial {\boldsymbol{u}}}}{{\partial t}} + \left({{\boldsymbol{u}} \cdot \nabla } \right){\boldsymbol{u}} = - \frac{1}{\rho }\nabla P + {\boldsymbol{g}} + c{\boldsymbol{u}}.$

The corrected PDEs of the system read

(30)

$\frac{{\partial {\boldsymbol{u}}}}{{\partial t}} + \nabla P = \frac{{\Delta t}}{2}\frac{{{\partial ^2}{\boldsymbol{u}}}}{{\partial {t^2}}} + c{\boldsymbol{u}}.$

Within potential theory, the combined free surface boundary condition is transformed into

(31)

$\frac{{{\partial ^2}\Phi }}{{\partial {t^2}}} + g\frac{{\partial \Phi }}{{\partial z}} - \frac{{\Delta t}}{2}\frac{{{\partial ^3}\Phi }}{{\partial {t^3}}} - c\frac{{\partial \Phi }}{{\partial t}} = 0\;\;\;\;\;\;{\rm{on}}\;\;z = 0.$

The dispersion relationship is thus rewritten as

(32)

$ - \frac{{{\rm i}\Delta t}}{2}{\omega ^3} - {\omega ^2} +{\rm i}c\omega + gk\tanh kh = 0.$

In Eq. (32), one can find that the dispersion relationship is identical to the one of a monochromatic linear wave when

$c = \Delta t{\omega ^2}/2$

. This indicates that by reasonably selecting the strength of the source term, the numerical dissipation effects can be greatly declined or even be nearly eliminated.

By letting

$k = {k_0} + \delta $

and adopting Eq. (26), Eq. (32) gives the solution as follows

(33)

$k = {k_0} + \frac{{\rm i}}{{2{C_{\rm g}}}}\left({ - c + \frac{{\Delta t}}{2}\omega _0^2} \right).$

The expression of velocity potential and surface elevation are eventually written as

(34a)

$\tag{34a}\Phi = - a\frac{{{\rm i}g}}{{{\omega _0}}}{{\rm e}^{\frac{1}{{2{C_{\rm g}}}}\left({c - \;\frac{1}{2}\rm\Delta t{\omega _0}^2} \right)x}}\frac{{\cosh {k_0}\left({z + h} \right)}}{{\cosh {k_0}h}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}},$

(34b)

$\tag{34b}\eta = a{{\rm e}^{\frac{1}{{2{C_{\rm g}}}}\left({c - \;\frac{1}{2}\rm\Delta t{\omega _0}^2} \right)x}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}}.$

It is worth noting that numerical dissipation is helpful to the stability of the solution. Therefore, the magnitude of the source term should not be too strong, otherwise unphysical energy boost will be induced and numerical instability occurs. In addition, the nonlinear and high-order terms of Eq. (9) can also affect the dissipation behavior of the corrected PDEs. Although the influence is not significant for small-amplitude waves, it cannot be neglected once the wave amplitude becomes large. For above reasons, Eq. (34) can be further written as

(35a)

$\tag{35a}\Phi = - a\frac{{{\rm i}g}}{{{\omega _0}}}{{\rm e}^{\left({\frac{c}{{2{C_g}}} - \varepsilon } \right)x}}\frac{{\cosh {k_0}\left({z + h} \right)}}{{\cosh {k_0}h}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}},$

(35b)

$\tag{35b}\eta = a{{\rm e}^{\left({\frac{c}{{2{C_g}}} - \varepsilon } \right)x}}{{\rm e}^{{\rm i}\left({{k_0}x - {\omega _0}t} \right)}}.$

where

$\varepsilon $

represents the contribution of numerical dissipation, which needs not to be

$\Delta t{\omega _0}^2/4{C_g}$

if nonlinear and high order effects are taken into account.

4 Numerical tests

Less

This section investigates numerical dissipation phenomena during wave propagation by generating a series of waves in NWTs. To verify the theory and method proposed in this work, unless otherwise specified, the NWTs utilize the same methods and parameters as those in Section 2. To be specific, the NWTs solve Euler equations with FDM and staggered grids, convective terms and the pressure gradient are discretized with CDS, the velocity-pressure system is de-coupled with a projection method, and time advancement is realized by a forward Euler scheme. It is worth noting that the instability induced by discretizing convective terms with CDS can be avoided by a typical deferred-correction approach (Ferziger and Peric, 2012). Wave generation is accomplished by directly assigning the boundary condition of the left-side boundary of the NWT with analytical solutions of corresponding waves. A conserved absorption method (Hu et al., 2015) is adopted at the right end of the tank for wave damping. The length of the wave absorption region equals the wavelength of the generated waves. A schematic sketch of the NWT is illustrated in Fig. 2.

4.1 Numerical dissipation of linear waves

In this test, small-amplitude linear waves are generated and propagating in a long NWT where the convective terms are discretized with CDS. The length and water depth of the NWT are 180 m and 0.5 m, respectively. Three groups of tests are conducted, with wave numbers (k) of 1.112 m^–1, 2.223 m^–1 and 3.335 m^–1, respectively. Wave amplitudes (a) of these tests range from 0.003 m to 0.015 m, and time advancing intervals change from 0.000 5 s to 0.005 s. The x-th and z-th grid number are 4 000 and 25, respectively.

Figure 3 shows contours of an example wave (k=2.223 m^–1, a=0.009 m) simulated with various time steps. Noted that the wave absorption region offers additional damping effect which is not of interest here, all contours within the wave absorption region are truncated and erased hereinafter. For comparison, the wave envelope shapes predicted by Eq. (28) are also displayed in red lines in Fig. 3. It can be clearly observed that the amplitudes of waves are reduced after long distance of travelling. The rate of reduction is relatively high in the vicinity of the wave generator, and becomes low when the wave is travelling away from the wave generator. Besides, with the growth of simulation time step, the damping effect becomes more and more distinct. The reduction curve of wave amplitudes can be well approximated by equation

$y = a \cdot \exp \left(- \;\dfrac{{\Delta t{\omega _0}^2}}{{4{C_g}}}x\right)$

which is extracted from Eq. (28).

To reveal the spatial-frequency variation of wave energy, a widely-used wavelet tool, i.e., the Gabor Transformation, is performed upon the wave contour results. The Gauss function is select as the window of Gabor Transformation. The Gabor transformation is formulated as (Daubechies, 1992)

(36)

$\left({Gf} \right)\left({b,k} \right) = \int {f\left(x \right)} {g^ * }\left({x - b} \right){{\rm e}^{ -{\rm i}kx}}{\rm d} x,$

where

$\,f\left(x \right)$

is the signal (wave contour) and

$ * $

denotes complex conjugate. Function

$g\left(x \right)$

is the window of Gabor transformation and is taken as

$1/\sqrt {4{\text π}{\text{α} _0}} \cdot \exp \left({ - {x^2}/4{\text{α} _0}} \right)$

in this article. The Gabor Transformations of the generated wave contours are shown in Fig. 4. The wave number (k=2.223 m^–1) of background carrier waves is also drawn with white dash-dot lines in Fig. 4. It is found that the wave energy shows clear reduction during propagation. With the increasing of the time advancing step, the speed of energy dissipation presents a remarkable boost. One can also find in Fig.4 that no frequency shift occurs after long distance of propagation, indicating that the key behavior of the numerical algorithm is numerical dissipation rather than dispersion.

Figure 5 displays numerical-predicted wave amplitudes at various positions, and contrasts the results against the wave envelope outlines calculated by Eq. (28). From Fig. 5, one can observe that Eq. (28) is capable of giving good approximation on the wave amplitude spatial distributions calculated with different time steps. However, it should also be noted that Eq. (28) slightly over-predicts the effect of numerical dissipation, thus under-estimating the wave amplitudes at the end of the tank. In addition, the deviation between numerical results and Eq. (28) becomes more evident with the growth of distance from the wave generator.

The propagation of waves with other amplitudes shows similar patterns as are illustrated by Fig. 3, Fig. 4 and Fig. 5, which are therefore not listed in the paper. By defining the wave amplitude reduction factor as

${R_W} = {a_1}/{a_0}$

, where

${a_0}$

and

${a_1}$

represents wave amplitudes nearby the wave generator and in front of the wave absorption region respectively. Figure 6 gives an overall description on the effect of numerical dissipation with various parameters. In Fig. 6, all

${R_W}$

obtained with various time steps, wave numbers and initial wave amplitudes (

${a_0}$

) are displayed together in a scatter diagram. One can find that

${R_W}$

is significantly impacted by length of time steps and wave numbers, whereas initial amplitudes barely affect

${R_{\rm W}}$

. The theoretical prediction by Eq. (28) is also plotted in Fig. 6, showing good agreement with the numerical simulation in a global sense although deviation exists.

As is proved by Eq. (35), the strength of the source term should be

$c = 2{C_{\rm g}}\varepsilon $

to completely eliminate the influence of numerical dissipation.

$\varepsilon $

can be approximated by

$\Delta t{\omega _0}^2/4{C_{\rm g}}$

only when nonlinear effects and high-order truncation error are ignored. In realistic applications, one can calculate

$\varepsilon $

in advance by simulating the free dissipation of waves with no compensation for numerical damping, which certainly accounts for nonlinear and high-order effects.

$c$

can subsequently be determined using formula

$c = 2{C_{\rm g}}\varepsilon $

. Following this idea, we adopt source terms with strengths listed in Table 1 to compensate for the numerical damping effect. Figure 7 displays the simulation results of the example wave (k=2.223 m^–1, a=0.009 m), where one can observe an apparent decrease of the numerical dissipation effect. In Fig. 8, the wave energy is barely changed after long term of travel, indicating that the influence of numerical damping is tiny. Figure 9 summarizes the wave amplitude reduction factors of all simulation results. All wave amplitude reduction factors calculated with source terms are close to 1, which demonstrates that the source terms works well for waves of various amplitudes, wave numbers and time step sizes.

4.2 Numerical dissipation of nonlinear Stokes waves

This subsection investigates the influence of numerical dissipation upon waves when nonlinear terms are taken into account. In tests of this subsection, 5-order Stokes waves are simulated under a water depth of 0.5 m. The wave number is 2.223 m^–1, the time step size is 0.003 2 s and wave heights are 0.06 m, 0.08 m and 0.1 m, respectively. Figure 10 shows the wave outlines of various moments when no source terms are used. The wave envelope predicted by Eq. (28) is also plotted in red dashed lines. Figure 10 reveals that numerical dissipation of nonlinear Stokes waves is stronger than that of small-amplitude linear waves.

The spatial energy distribution of the Stokes waves is displayed in Fig. 11. Similar to Fig. 4, clear energy consumption can be observed in Fig. 11. Besides, one can find deviation between actual wave frequencies (or wave numbers) and the intended ones (i.e., white dashed lines). With the growth of distance from wave generators, the deviation seems more and more palpable. The wavelengths of the simulated waves seem to slowly increase during propagation, which is more probably induced by a numerical dispersion effect. This indicates that numerical dispersion may also have to be considered in addition to numerical dissipation when simulating large-amplitude nonlinear waves.

According to Eq. (35), the strengths of source terms are taken as c=0.055, c=0.058 and c=0.063 with respect to Fig. 10a, Fig. 10b and Fig. 10c, respectively. Taking into account the influence of source terms, the simulated wave contours are drawn in Figure 12. Figure 13 shows the spatial energy of the generated waves. Figure 12 and Fig. 13 reveal that the existence of source terms can effectively balance the numerical dissipation during nonlinear wave simulation.

During wave simulation, vertical grid dimension is a crucial reason affecting the quality of the generated waves. Figure 14 simulates the same 5-order Stokes wave (k=2.223 m^–1 and wave height=0.06 m) with various vertical grid dimensions and no compensating source term. Figure 15 simulates the same waves as Fig. 14, but with compensating source terms equaling 0.055, 0.028 and 0.027 for subplot (a), (b) and (c), respectively. Table 2 lists the wave heights of Fig. 14 and Fig. 15 at position x=5 m and x=10 m, and contrast the simulation results against theory. Figure 14 and Table 2 reveal that denser vertical grid size can improve the simulation precision, yet the rate of improvement becomes slow once the vertical grid size has been enough. Figure 15 and Table 2 show that the compensating source term works well under various vertical grid densities.

In order to test how the proposed method works for waves with very large steepness, we try to generate a high 5-order Stokes wave whose height almost reaches the breaking limit. The wave number and wave height are 2.223 m^–1 and 0.32 m, respectively. Fig. 16 shows the wave outlines at various moments computed a without any source term and Fig.16b with a source term whose strength equals 0.092. Table 3 lists wave heights at various locations and their relative error with respect to the target value. One can find that large-height wave are more easily affected by numerical dissipation, and the reduction of wave height is significant. It should be mentioned that the wave-height reduction rate can be declined by using various techniques, to name but a few, using denser grids or high-order difference scheme. However, these methods are not of concern in this paper, and we reserve the large dissipation effect to give a clearer observation on the dissipation reduction ability of the source term. By using the compensating source term, the dissipation effect is significantly weakened which can be seen from Fig. 16b and Table 3, despite that an error of about 20% still exists. One can try larger source strength to make the result closer to the theory. However, this is not recommended since larger source strength brings more instability factors to the solution, and the solution may diverge easily under certain situations.

4.3 Numerical dissipation of irregular waves

Irregular waves consist of wave components with various frequencies, wave numbers and amplitudes. The numerical dissipation of irregular waves is more complex due to the fact that the strength of damping is different for each wave component. Moreover, the focusing of wave components may induce a huge wave crest, and therefore lead to significant enhancement in nonlinearity. Without loss of generality, this subsection studies the numerical dissipation behavior of irregular waves by simulating the P-M spectrums in a 90 m-long numerical wave flume with a water depth of 0.5 m. The zero-crossing period of the spectrum is

${T_0} = $

1.5 s, and the significant wave heights are

${H_S}$

= 0.01 m, 0.03 m and 0.05 m, respectively.

Figure 17 displays surface elevations of the irregular wave with

${H_S}$

=0.05 m at x=0 m (i.e., at the wave generating boundary), x=9 m and x=18 m. From Fig. 17, one can clearly observe the energy damping phenomenon during wave propagation. Furthermore, the damping of large-amplitude components is severer. For instance, the maximum wave height of subplot (c) is much smaller than that of subplot (a), indicating that nonlinearity seems to strengthen the numerical damping effect. The conclusion is also proved by Fig. 18, where energy spectrums of the simulated waves at various locations are plotted. In Fig. 18, the simulated spectrums show well agreement with the target spectrums nearby the wave generator. However, the entire energy of the wave spectrum is rapidly declined with the growth of x coordinates. The rate of decline becomes faster under larger

${H_S}$

values. It is worthy of noting that the simulated wave spectrums become narrower than the target spectrums due to the reason that high-frequency wave components are more sensitive to numerical dissipation, thus making the high-frequency portion of the spectrum decline significantly.

In comparison with regular waves, the compensation for numerical dissipation of irregular waves requires additional treatment. The relationship between the energy of the simulated waves and the simulated spectrums reads

(37)

$E' = \int_0^\infty {S'(\omega){\rm d}\omega } = \sum\limits_{i = 1}^N {\frac{1}{2}A_i^2},$

where

$E'$

and

$S'$

are wave energy and wave spectrum with the influence of source terms, respectively,

${A_i}$

is the amplitude of the i-th wave component and N the number of composing waves. Using Eq. (35), the right side of Eq. (37) is written as

(38)

$\sum\limits_{i = 1}^N {\frac{1}{2}A_i^2} = \sum\limits_{i = 1}^N {\frac{1}{2}{a_i}^2{{\rm e}^{\left({\frac{c}{{{C_{gi}}}} - 2{\varepsilon _i}} \right)x}}} .$

Under the hypothesis of narrow-band spectrum, Eq. (38) is further simplified as

(39)

$\sum\limits_{i = 1}^N {\frac{1}{2}A_i^2} = {{\rm e}^{\frac{c}{{{C_g}}}x}}\sum\limits_{i = 1}^N {\frac{1}{2}{a_i}^2{{\rm e}^{ - 2{\varepsilon _i}x}}} = {{\rm e}^{\frac{c}{{{C_g}}}x}}\int_0^\infty {S(\omega)d\omega } = {{\rm e}^{\frac{c}{{{C_g}}}x}}E,$

where

$E$

and

$S$

are wave energy and wave spectrum simulated without using any source term (only affected by numerical dissipation). The main principle for the determination of the strength of the source term is to ensure the energy of simulated waves to be equal to the energy of target spectrums, namely

$E' = {E_0}$

. Therefore, by combining Eq. (37) and Eq. (39), the source strength if finally decided by

(40)

$c = \frac{{{C_g}}}{x}\ln \frac{{{E_0}}}{E} = \frac{{{C_g}}}{x}\ln \frac{{\displaystyle\int_0^\infty {{S_0}(\omega){\rm d}\omega } }}{{\displaystyle\int_0^\infty {S(\omega){\rm d}\omega } }},$

where

${S_0}$

and

${E_0}$

are the target wave spectrum and its energy.

By applying Eq. (40), the strengths of compensating source terms are c=0.011 2, c=0.017 8 and c=0.057 2 for subplots Fig. 18a, Fig. 18b and Fig. 18c, respectively. Appending these source terms to Euler equations as is expressed in Eq. (29), the irregular waves of this subsection are regenerated, and the spectrums of various waves at different locations are depicted in Fig. 19. In Fig. 19, the spectrums at various positions exhibit similar shapes. By comparing Fig. 19 against Fig. 18, it is clear that numerical dissipation is effectively weakened.

5 Conclusions

Less

On the basis of above simulations and discussions, the following conclusions can be drawn.

(1) Equation (28) can be used to describe the dissipation law of small-amplitude linear waves. Yet palpable deviation can be found if Eq. (28) is utilized for description of the dissipation of nonlinear Stokes waves.

(2) Compared with long waves, short waves are more likely to be affected by numerical dissipation effects, namely are easier to be dissipated.

(3) Due to Conclusion (2), the spectrum shapes are changed during wave propagation.

(4) Due to Conclusion (2), large-amplitude nonlinear Stokes waves are easier to be damped in contrast with small-amplitude linear waves, as high-frequency terms (i.e., nonlinear terms) are sensitive to numerical dissipation.

(5) The proposed method is capable of effectively reducing the effect of numerical dissipation during the simulation of small-amplitude linear waves, nonlinear Stokes wave and irregular wave spectrums.

Funding

Less

The National Natural Science Foundation of China under contract No. 51609101 and 51909103; the Natural Science Foundation of Fujian Province of China under contract Nos 2017J01701, 2017J05085 and 2018J05090; the Outstanding Young University Scientific Research Talents Cultivation Plan of Fujian Province of China.

References

Less

Abbasnia A, Ghiasi M. 2015. Fully nonlinear wave interaction with an array of truncated barriers in three dimensional numerical wave tank. Engineering Analysis with Boundary Elements, 58: 79–85, doi: 10.1016/j.enganabound.2015.03.015

Abbasnia A, Ghiasi M, Abbasnia A. 2017. Irregular wave transmission on bottom bumps using fully nonlinear NURBS numerical wave tank. Engineering Analysis with Boundary Elements, 82: 130–140

Abbasnia A, Soares C G. 2018. Transient fully nonlinear ship waves using a three-dimensional NURBS numerical towing tank. Engineering Analysis with Boundary Elements, 91: 44–49, doi: 10.1016/j.enganabound.2018.03.011

Alvarado-Rodríguez C E, Klapp J, Sigalotti L D G, et al. 2017. Nonreflecting outlet boundary conditions for incompressible flows using SPH. Computers & Fluids, 159: 177–188

Anbarsooz M, Passandideh-Fard M, Moghiman M. 2013. Fully nonlinear viscous wave generation in numerical wave tanks. Ocean Engineering, 59: 73–85, doi: 10.1016/j.oceaneng.2012.11.011

Beljadid A, LeFloch P G, Mishra S, et al. 2017. Schemes with well-controlled dissipation. hyperbolic systems in nonconservative form. Communications in Computational Physics, 21(4): 913–946

Bihs H, Kamath A, Chella M A, et al. 2016. A new level set numerical wave tank with improved density interpolation for complex wave hydrodynamics. Computers & Fluids, 140: 191–208

Cao Zhiwei, Liu Zhifeng, Wang Xiaohong, et al. 2017. A dissipation-free numerical method to solve one-dimensional hyperbolic flow equations. International Journal for Numerical Methods in Fluids, 85(4): 247–263, doi: 10.1002/fld.4383

Daubechies I. 1992. Ten Lectures on Wavelets. Philadelphia, PA: Society for Industrial and Applied Mathematics

De Paulo G S, Tomé M F, McKee S. 2007. A marker-and-cell approach to viscoelastic free surface flows using the PTT model. Journal of Non-Newtonian Fluid Mechanics, 147(3): 149–174, doi: 10.1016/j.jnnfm.2007.08.003

Dean R G, Dalrymple R A. 1991. Water Wave Mechanics for Engineers & Scientists (Vol. 2). Singapore: World Scientific Publishing Company

Elhanafi A, Macfarlane G, Fleming A, et al. 2017. Experimental and numerical measurements of wave forces on a 3D offshore stationary OWC wave energy converter. Ocean Engineering, 144: 98–117, doi: 10.1016/j.oceaneng.2017.08.040

Ferziger J H, Peric M. 2012. Computational Methods for Fluid Dynamics. Berlin Heidelberg: Springer

Hasan S A, Sriram V, Selvam R P. 2018. Numerical modelling of wind-modified focused waves in a numerical wave tank. Ocean Engineering, 160: 276–300, doi: 10.1016/j.oceaneng.2018.04.044

Hu Zhe, Tang Wenyong, Xue Hongxiang, et al. 2015. Numerical simulations using conserved wave absorption applied to Navier–Stokes equation model. Coastal Engineering, 99: 15–25, doi: 10.1016/j.coastaleng.2015.02.007

Hu Zhe, Tang Wenyong, Xue Hongxiang, et al. 2017. Numerical study of rogue wave overtopping with a fully-coupled fluid-structure interaction model. Ocean Engineering, 137: 48–58, doi: 10.1016/j.oceaneng.2017.03.022

Li Zhao, Zhang Yufei, Chen Haixin. 2015. A low dissipation numerical scheme for implicit large eddy simulation. Computers & Fluids, 117: 233–246

Liu Xin, Lin Pengzhi, Shao Songdong. 2015. ISPH wave simulation by using an internal wave maker. Coastal Engineering, 95: 160–170, doi: 10.1016/j.coastaleng.2014.10.007

Ma Z H, Causon D M, Qian L, et al. 2016. Numerical investigation of air enclosed wave impacts in a depressurised tank. Ocean Engineering, 123: 15–27, doi: 10.1016/j.oceaneng.2016.06.044

Nazari F, Mohammadian A, Charron M. 2015. High-order low-dissipation low-dispersion diagonally implicit Runge–Kutta schemes. Journal of Computational Physics, 286: 38–48, doi: 10.1016/j.jcp.2015.01.020

Panicker P G, Goel A, Iyer H R. 2015. Numerical modeling of advancing wave front in dam break problem by incompressible navier-stokes solver. Aquatic Procedia, 4: 861–867, doi: 10.1016/j.aqpro.2015.02.108

Park J C, Uno Y, Sato T, et al. 2004. Numerical reproduction of fully nonlinear multi-directional waves by a viscous 3D numerical wave tank. Ocean Engineering, 31(11–12): 1549–1565

Saincher S, Banerjeea J. 2015. Design of a numerical wave tank and wave flume for low steepness waves in deep and intermediate water. Procedia Engineering, 116: 221–228, doi: 10.1016/j.proeng.2015.08.394

Schillaci E, Jofre L, Balcázar N, et al. 2016. A level-set aided single-phase model for the numerical simulation of free-surface flow on unstructured meshes. Computers & Fluids, 140: 97–110

Schranner F S, Domaradzki J A, Hickel S, et al. 2015. Assessing the numerical dissipation rate and viscosity in numerical simulations of fluid flows. Computers & Fluids, 114: 84–97

Soares D Jr. 2019. A simple explicit-implicit time-marching technique for wave propagation analysis. International Journal of Computational Methods, 16(1): 1850082, doi: 10.1142/S0219876218500822

Appendix

Less

Year 2020 volume 39 Issue 1

PDF

Cite this Article

BibTeX

Article Info

doi: 10.1007/s13131-019-1524-1

Receive Date：2018-08-29
Online Date：2026-03-27
Published：2020-01-25

Article Data

Affiliations

History

Received：2018-08-29
Accepted：2019-01-06

Funding

The National Natural Science Foundation of China under contract No. 51609101 and 51909103; the Natural Science Foundation of Fujian Province of China under contract Nos 2017J01701, 2017J05085 and 2018J05090; the Outstanding Young University Scientific Research Talents Cultivation Plan of Fujian Province of China.

Affiliations

¹ Key Laboratory of Ships and Ocean Engineering of Fujian Province, School of Marine Engineering, Jimei University, Xiamen 361021, China

² Deep Sea Technology Research Center, School of Engineering, West Lake University, Hangzhou 310024, China

³ Hadal Science and Technology Research Center, College of Marine Sciences, Shanghai Ocean University, Shanghai 201306, China

Corresponding:

* E-mail: zhangxy@jmu.edu.cn

References

Share

https://castjournals.cast.org.cn/joweb/aos/EN/10.1007/s13131-019-1524-1

Share to

Scan QR to access full text

Cite this article

BibTeX

Citations

表12种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

Table 1. Strength of source terms, used for suppressing the numerical dissipation of small-amplitude linear waves

Time step size/s	Source term strength c/s^–1
Time step size/s	k=1.112 m^–1	k=2.223 m^–1	k=3.335 m^–1
0.000 5	0.000 84	0.002 66	0.004 62
0.001 0	0.002 22	0.007 06	0.012 26
0.002 5	0.005 92	0.018 87	0.032 76
0.005 0	0.012 39	0.039 45	0.068 50

Table 2. Wave heights of 5-order Stokes waves (k=2.223 m^-1 and target wave height=0.06 m), simulated with various vertical grid numbers, with and without the compensating source term

Vertical grid number	Simulation without compensating source term				Simulation with compensating source term
	Absolute wave height/m		Error against theory/%		Absolute wave height/m		Error against theory/%
	at x=5 m	at x=10 m	at x=5 m	at x=10 m	at x=5 m	at x=10 m	at x=5 m	at x=10 m
25	0.051 4	0.046 7	14.3	22.2	0.055 4	0.053 3	7.7	11.2
50	0.054 3	0.052 7	9.5	12.2	0.056 6	0.056 4	5.7	6.0
75	0.055	0.052 9	8.3	11.8	0.056 9	0.056 7	5.2	5.5

Table 3. Wave heights of 5-order Stokes waves (k=2.223 m^–1 and target wave height=0.32 m) at various positions, simulated with and without the compensating source term

Distance from wave generator/m	No source term		With compensating source term
Distance from wave generator/m	Wave height/m	Error against theory/%	Wave height/m	Error against theory/%
0.1	0.308	3.8	0.311	2.8
5	0.173	45.9	0.256	20.0
10	0.123	61.6	0.250	21.2
15	0.095	70.3	0.253	20.9

Fig. 1. Sketch of the staggered grid.

Fig. 2. Sketch of the numerical wave tank.

Fig. 3. Contours of small-amplitude linear waves (k=2.223 m^–1, a=0.009 m), simulated with time steps equaling (a) 0.000 5 s, (b) 0.001 s, (c) 0.002 5 s and (d) 0.005 s.

Fig. 4. Gabor transformation of contours of small-amplitude linear waves (k=2.223 m^–1 in white dash-dot line, a=0.009 m), simulated with time steps of 0.000 5 s (a) , 0.001 s (b), 0.002 5 s (c) and 0.005 s (d).

Fig. 5. Spatial wave-amplitude distributions of small-amplitude linear waves (k=2.223 m^–1, a=0.009 m), predicted by numerical simulation (lines) and theoretical formula Eq. (28) (marker points).

Fig. 6. Wave amplitude reduction factors of small-amplitude linear waves, simulated with various amplitudes, wave numbers and time steps.

Fig. 7. Contours of small-amplitude linear waves (k=2.223 m^–1, a=0.009 m), simulated with time steps equaling 0.000 5 s (a), 0.001 s (b), 0.002 5 s (c) and 0.005 s (d). Black solid lines denote wave contours simulated with source terms, and red filled circles wave amplitudes simulated without any source term.

Fig. 8. Gabor transformation of contours of small-amplitude linear waves (k=2.223 m^–1 white dash-dot line, a=0.009 m), simulated considering the influence of source terms, with time steps of 0.000 5 s (a), 0.001 s (b), 0.002 5 s (c) and 0.005 s (d).

Fig. 9. Wave amplitude reduction factors of small-amplitude linear waves, simulated with (scatter points) and without (lines) the influence of source terms.

Fig. 10. Wave outlines of 5-order Stokes waves with k=2.223 m^–1 and wave heights equal to 0.06 m (a), 0.08 m (b) and 0.10 m (c) . Red dashed lines denote envelopes calculated by Eq. (28).

Fig. 11. Gabor transformation of outlines of 5-order Stokes waves (k=2.223 m^–1 white dash-dot line), with wave heights of 0.06 m (a), 0.08 m (b) and 0.1 m (c).

Fig. 12. Wave outlines of 5-order Stokes waves with k=2.223 m^–1 and wave heights equal to 0.06 m (a), 0.08 m (b) and 0.10 m (c). Source terms are involved during simulation.

Fig. 13. Gabor transformation of outlines of 5-order Stokes waves (k=2.223 m^–1 white dash-dot line), with wave heights of 0.06 m (a), 0.08 m (b) and 0.1 m (c). Source terms are involved during simulation.

Fig. 14. Wave outlines of 5-order Stokes waves with k=2.223 m^–1 and wave heights=0.06 m, with z-th grid number equaling 25 (a), 50 (b) and 75 (c).

Fig. 15. Wave outlines of 5-order Stokes waves with k=2.223 m^–1 and wave height=0.06 m, with z-th grid number equaling 25 (a), 50 (b) and 75 (c). Source terms are involved during simulation.

Fig. 16. Wave outlines of 5-order Stokes waves with k=2.223 m^–1 and wave height=0.32 m, without (a) and with (b) the compensating source term.

Fig. 17. Surface histories of an irregular wave based on the P-M spectrum (${H_{\rm S}}$=0.05 m and ${T_0}$=1.5 s), at x=0 m (a), x=9 m (b) and x=18 m (c) .

Fig. 18. Energy spectrums of irregular waves (${T_0}$=1.5 s) based on the P-M spectrum, with ${H_{\rm S}}$=0.01 m (a), ${H_{\rm S}}$=0.03 m (b) and ${H_{\rm S}}$=0.05 m (c).

Fig. 19. Energy spectrums of irregular waves (${T_0}$=1.5 s) based on the P-M spectrum involving the influence of source terms, with ${H_{\rm S}}$=0.01 m (a), ${H_{\rm S}}$=0.03 m (b) and ${H_{\rm S}}$=0.05 m (c).

Articles: Latest Articles; Most Read; Collections

Updates: Events; News; Multimedia

About: About Us

Contact

No. 86 Xueyuan South Road, Haidian District, Beijing

100081

010-62199257

qkjq@cast.org.cn

Copyright © 2025 China Association for Science and Technology. All rights reserved. For all open access content, the relevant licensing terms apply.
Sponsored by the Office of the Leading Group for Cybersecurity and Informatization of CAST, and supported by Science and Technology Review Publishing House