Accurate fast computation of steady two-dimensional surface gravity waves in arbitrary depth

TL;DR

This paper introduces an efficient, open-source algorithm for accurately computing steady two-dimensional surface gravity waves in any depth, suitable for various wave types and steepnesses, using conformal mapping and spectral methods.

Contribution

The paper presents a novel, fast algorithm with O(N log N) complexity for arbitrary depth wave computation, combining conformal mapping, Babenko equation, and Petviashvili's iterations.

Findings

Algorithm achieves high precision for Stokes, cnoidal, and solitary waves.

Efficient for large steepnesses and arbitrary depths.

Open-source code enables easy verification and adaptation.

Abstract

This paper describes an efficient algorithm for computing steady two-dimensional surface gravity wave in irrotational motion. The algorithm complexity is O(N log N), N being the number of Fourier modes. The algorithm allows the arbitrary precision computation of waves in arbitrary depth, i.e., it works efficiently for Stokes, cnoidal and solitary waves, even for quite large steepnesses. The method is based on conformal mapping, Babenko equation rewritten in a suitable way, pseudo-spectral method and Petviashvili's iterations. The efficiency of the algorithm is illustrated via some relevant numerical examples. The code is open source, so interested readers can easily check the claims, use and modify the algorithm.

Tables

Table 1. Table 1. KdV parameter m 𝑚 m for H / d = 0.1 𝐻 𝑑 0.1 H/d\ =\ 0.1\, .

$L/d$	$100$	$\mathrm{1\hspace{0.17em}000}$	$\mathrm{10\hspace{0.17em}000}$	$\mathrm{100\hspace{0.17em}000}$	$\mathrm{1\hspace{0.17em}000\hspace{0.17em}000}$
$1-m$	$2.04\times {10}^{-11}$	$1.85\times {10}^{-118}$	$6.91\times {10}^{-1189}$	$3.63\times {10}^{-11893}$	$5.78\times {10}^{-118936}$

Table 2. Table 2. KdV parameter m 𝑚 m for H / d = 0.5 𝐻 𝑑 0.5 H/d\ =\ 0.5\, .

$L/d$	$50$	$60$	$70$	$80$	$90$	$100$
$1-m$	$8.06\times {10}^{-13}$	$1.77\times {10}^{-15}$	$3.87\times {10}^{-18}$	$8.47\times {10}^{-21}$	$1.86\times {10}^{-23}$	$4.06\times {10}^{-26}$

Table 3. Table 3. Various parameters for the highest computable waves.

$kd$	$\hat{\epsilon }$	${\epsilon }_{\text{b}}$	${\epsilon }_{\text{c}}$	$A\left({\epsilon }_{\text{c}}\right)$	$\vartheta \left({\epsilon }_{\text{c}}\right)$	${\theta }_{\text{max}}\left({\epsilon }_{\text{c}}\right)$
$\infty$	$0.443164$	$0.435907$	$0.4400$	$2.4366$	$2.4710$	${29.831}^{\circ }$
$2.0$	$0.422293$	$0.415166$	$0.4199$	$2.5651$	$2.5865$	${29.942}^{\circ }$
$1.5$	$0.389554$	$0.382626$	$0.3876$	$2.6336$	$2.6475$	${29.979}^{\circ }$
$1.0$	$0.315872$	$0.309415$	$0.3146$	$2.7535$	$2.7573$	${30.042}^{\circ }$
$0.9$	$0.294147$	$0.287847$	$0.2930$	$2.7714$	$2.7734$	${30.050}^{\circ }$
$0.8$	$0.269865$	$0.263760$	$0.2688$	$2.7656$	$2.7674$	${30.041}^{\circ }$
$0.7$	$0.243084$	$0.237227$	$0.2421$	$2.7525$	$2.7547$	${30.028}^{\circ }$
$0.6$	$0.213965$	$0.208426$	$0.2130$	$2.6964$	$2.7007$	${29.981}^{\circ }$
$0.5$	$0.182750$	$0.177626$	$0.1818$	$2.6216$	$2.6297$	${29.910}^{\circ }$
$0.4$	$0.149709$	$0.145132$	$0.1488$	$2.5396$	$2.5520$	${29.822}^{\circ }$
$0.3$	$0.115019$	$0.111166$	$0.1142$	$2.4542$	$2.4723$	${29.714}^{\circ }$
$0.2$	$0.078662$	$0.075763$	$0.0780$	$2.3636$	$2.3886$	${29.579}^{\circ }$
$0.106814$	$0.043100$	$0.041362$	$0.0427$	$2.3099$	$2.3399$	${29.488}^{\circ }$

Equations

$$⟨\eta ⟩\stackrel{\mathrm{def}}{:=}\frac{k}{2\pi }{\int }_{-\pi /k}^{\pi /k}\eta (x)\mathrm{d}\text{/}x=0.$$

$$\displaystyle\phi_{\,x\,x}\ +\ \phi_{\,y\,y}\$$

$$\displaystyle\phi_{\,y}\$$

$$\displaystyle\phi_{\,y}\ -\ \eta_{\,x}\,\phi_{\,x}\$$

$$\displaystyle 2\,g\,\eta\ +\ \phi_{\,x}^{\,2}\ +\ \phi_{\,y}^{\,2}\$$

$$2p+2gy+u^2+v^2=B,$$

$$B=⟨{\tilde{u}}^2+{\tilde{v}}^2⟩.$$

$$B=⟨{\bar{u}}^2⟩,$$

$$c_{\text{textsc}s}\stackrel{\mathrm{def}}{:=}-⟨\frac{1}{0pt}{\int }_{-0pt}^{\eta }u(x,y)\mathrm{d}\text{/}y⟩=\frac{\bar{\psi }-\tilde{\psi }}{0pt}.$$

$$c_{\text{textsc}e}\stackrel{\mathrm{def}}{:=}-⟨\bar{u}⟩=-⟨u(x,y=-d)⟩.$$

$$0⩽x⩽2\pi /k,-0pt⩽y⩽\eta (x),$$

$$0⩽\alpha ⩽2\pi c_{\text{textsc}e}/k{c}_{\text{textsc}r},-0pt{c}_{\text{textsc}s}/c_{\text{textsc}r}⩽\beta ⩽0,$$

$$\text{mathchar}22k\stackrel{\mathrm{def}}{:=}{c}_{\text{textsc}r}{c}_{\text{textsc}e}^{-1}k,\text{mathchar}22d\text{/}\stackrel{\mathrm{def}}{:=}{c}_{\text{textsc}s}{c}_{\text{textsc}r}^{-1}d,$$

$$\begin{gathered} \frac{w}{c_{\text{textsc}r}}=-{\left(\frac{\mathrm{d}z}{\mathrm{d}\text{/}\zeta }\right)}^{-1},\frac{u}{c_{\text{textsc}r}}=\frac{-x_{\alpha }}{x_{\alpha }^2+y_{\alpha }^2}, \\ \frac{v}{c_{\text{textsc}r}}=\frac{-y_{\alpha }}{x_{\alpha }^2+y_{\alpha }^2},\frac{u^2+v^2}{c_{\text{textsc}r}^2}=\frac{1}{x_{\alpha }^2+y_{\alpha }^2} \end{gathered}$$

$$x_{\alpha }=y_{\beta }=\frac{-c_{\text{textsc}r}u}{u^2+v^2},y_{\alpha }=-x_{\beta }=\frac{-c_{\text{textsc}r}v}{u^2+v^2}.$$

$$\frac{B}{c_{\text{textsc}r}^2}=\frac{k}{2\pi }{\int }_{-\frac{\pi c_{\text{textsc}e}}{k{c}_{\text{textsc}r}}}^{\frac{\pi c_{\text{textsc}e}}{k{c}_{\text{textsc}r}}}\frac{x_{\alpha }}{x_{\alpha }^2+y_{\alpha }^2}\mathrm{d}\text{/}\alpha =\frac{c_{\text{textsc}e}}{c_{\text{textsc}r}}\frac{\text{mathchar}22k}{2\pi }{\int }_{-\pi /\text{mathchar}22k}^{\pi /\text{mathchar}22k}\frac{x_{\alpha }}{x_{\alpha }^2+y_{\alpha }^2}\mathrm{d}\text{/}\alpha ,$$

$$\left\llangle(\cdots)\right\rrangle\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \frac{\mathchar 22\relax\mkern-9.0muk}{2\,{\pi}}\;\int_{-{\pi}/\mathchar 22\relax\mkern-9.0muk}^{{\pi}/\mathchar 22\relax\mkern-9.0muk}\,(\cdots)\ \mathrm{d}\/\alpha\,,$$

$$\displaystyle\left\llangle\widetilde{(\cdots)}\right\rrangle\$$

$$\displaystyle\left\llangle\overline{(\cdots)}\right\rrangle\$$

$$\displaystyle\left<\widetilde{(\cdots)}\right>\$$

$$⟨\overline{(\cdots )}⟩$$

$$x=c_{\text{textsc}r}{c}_{\text{textsc}e}^{-1}\alpha +X(\alpha ,\beta ),y=c_{\text{textsc}r}{c}_{\text{textsc}e}^{-1}(\text{/}\beta +\text{mathchar}22d\text{/}\text{/})-0pt+Y(\alpha ,\beta ),$$

$$\bar{Y}(\alpha )\stackrel{\mathrm{def}}{:=}Y(\alpha ,-\text{mathchar}22d\text{/})=0,\tilde{Y}(\alpha )\stackrel{\mathrm{def}}{:=}Y(\alpha ,0)=\tilde{y}+0pt(1-c_{\text{textsc}s}/c_{\text{textsc}e}).$$

$$X(0,\beta )=0,X(2\pi /\text{mathchar}22k,\beta )=0,$$

$$X(\alpha ,\beta )=$$

$$=$$

$$Y(\alpha ,\beta )=$$

$$=$$

$$\tilde{X}(\alpha )=\cos \left[\text{/}\text{mathchar}22d\text{/}{\partial }_{\alpha }\text{/}\right]\bar{X}(\alpha ),$$

$$\bar{X}(\alpha )=\sec \left[\text{/}\text{mathchar}22d\text{/}{\partial }_{\alpha }\text{/}\right]\tilde{X}(\alpha ).$$

$$\tilde{Y}(\alpha )=\mathcal{H}\left\{\tilde{X}(\alpha )\right\},\mathcal{H}\stackrel{\mathrm{def}}{:=}\tan \left[\text{mathchar}22d\text{/}{\partial }_{\alpha }\right],$$

$$\mathcal{H}\left\{{\mathrm{e}}^{\mathrm{i}\kappa \alpha }\right\}=\text{mathchar}22d\text{/}\frac{\partial {\mathrm{e}}^{\mathrm{i}\kappa \alpha }}{\partial \alpha }+\frac{\text{mathchar}22d{\text{/}}^3}{3}\frac{{\partial }^3{\mathrm{e}}^{\mathrm{i}\kappa \alpha }}{\partial {\alpha }^3}+\frac{2\text{mathchar}22d{\text{/}}^5}{15}\frac{{\partial }^5{\mathrm{e}}^{\mathrm{i}\kappa \alpha }}{\partial {\alpha }^5}+\cdots =\mathrm{i}\tanh (\text{/}\kappa \text{/}\text{mathchar}22d\text{/}){\mathrm{e}}^{\mathrm{i}\kappa \alpha }.$$

$$\left\llangle\,\bar{x}\,\right\rrangle\ =\ c_{\,\textsc{r}}\,c_{\,\textsc{e}}^{\,-1}\left\llangle\,\alpha\,\right\rrangle\ +\ \left\llangle\,\bar{X}\,\right\rrangle\ =\ {\pi}\,k^{\,-1}\ +\ \left\llangle\,\bar{X}\,\right\rrangle\,.$$

$$\left\llangle\,\bar{x}\,\right\rrangle\ =\ \left<\,-\,c_{\,\textsc{e}}^{\,-1}\,x\,\bar{\phi}_{\,x}\,\right>\,=\ {\pi}\,k^{\,-1}\,.$$

Full text

Key words and phrases: Periodic waves; cnoidal waves; Stokes wave; Babenko equation; Petviashvili scheme; spectral methods

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Definitions and notations

3. 3 Conformal mapping

4. 3.1 Conformal averaging operator

5. 3.2 Resolution of the conformal mapping

6. 3.3 Averaging of dependent functions

7. 3.4 Mean level condition

8. 3.5 Celerities

9. 4 Babenko equations

10. 5 Numerical resolution

11. 5.1 User provided parameters

12. 5.2 Computational parameters

13. 5.3 Classical Petviashvili method (CPM)

14. 5.4 Computation of the unknown parameters

15. 5.5 Post processing

16. 6 Numerical examples

17. 6.1 Deep water

18. 6.2 Finite depth

19. 6.3 Shallow water

20. 6.4 Comparison with KdV cnoidal wave

21. 6.5 Solitary waves

22. 7 Remarks on the highest computable waves

23. 8 Discussion

24. Acknowledgments

25. A Potential and stream function in a ‘fix’ frame of reference

26. B Integral quantities

27. C Relations between averaged quantities in physical and conformal planes

28. D Velocity and pressure fields in the fluid

Introduction

Many physical phenomena and mathematical problems related to surface gravity water waves remain unknown, not well-understood or unsolved, even in the ‘simple’ case of traveling waves of permanent form in two-dimensional irrotational motion [10, 18]. Traveling waves are of special interest because complex sea states are often described as superposition and interaction of such waves. Surveys of analytical and numerical models, their limitations and open questions can be found in dedicated articles [20, 24, 25, 29, 44] and books [17, 38, 46].

Since exact analytic solutions for irrotational steady surface gravity waves are still unknown, and likely will never be known, only analytic or numerical approximations are accessible. Simple analytic approximations are interesting for physical insights, but they are of limited accuracy. Even formal analytic solutions in terms of small parameter expansions have limited accuracy because they generally converge slowly [41], when they converge (shallow water expansions are divergent for all amplitudes [28]). Thus, their numerical calculation suffer large truncation errors, and are prone to important accumulation of round-off and cancelation errors. Even when a simple analytic approximation would be sufficient for a given application, its computation may be practically intractable. An example is KdV cnoidal wave analytic solution that cannot be easily computed for very long waves, as demonstrated in the present paper. Therefore, only numerical approximations of the original equations can provide highly accurate solutions that are necessary for many applications. For instance, for stability analysis and interactions of traveling waves using accurate numerical models, a too crude (about six digits accuracy, say) initial condition for traveling wave may lead to incorrect behaviours, specially for long-time simulations. A lack of initial accuracy may then lead to erroneous physical interpretations. Another example is the numerical investigation of mathematical conjectures. With arbitrary precision computations, one can check if a conjecture is likely true or not, or can formulate new conjectures worthy investigations. Indeed, some questions of mathematical interest (e.g., rigorous proofs of unicity and monotonicity regarding velocity, acceleration and pressure fields) remain open for all waves and not only for very large amplitudes (see the papers in [18] for reviews).

Several algorithms have been proposed in the literature for the computation of steady surface waves solutions of the irrotational Euler equations [38]. The focus was more on the computation of the almost highest gravity waves (see, e.g., [7, 35, 47]) or ‘exotic’ capillary-gravity waves (e.g., [16, 46]), than the computation of arbitrary wavelengths. Actually, none of these algorithms are capable of computing long waves in shallow water, not even the ones of small amplitude, because they are too demanding. Indeed, all these algorithms lead to the resolution of a discrete system of nonlinear equations, this system being large for long waves. The resolution of this system is generally performed with (the like of) Newton or, better, Levenberg–Marquardt iterations [34]. Though robust and effective, these methods are computationally very demanding because each iteration requires $\mathcal{O}\,(N^{\,3})$ operations, $N$ being the number of unknowns (e.g. the Fourier coefficient for periodic water waves). Indeed, though their theoretical minimum complexity is $\mathcal{O}\,(N^{\,2})\,$, usually Levenberg–Marquardt iterations use direct solvers that have $\mathcal{O}(N^{\,3})$ complexity. In the simplest case, they rely on Cholesky decomposition which needs $\mathcal{O}\,(N^{\,3})$ operations. In most numerically robust version of the Levenberg–Marquardt method, its implementations relies on rank-revealing QR algorithm, which is also $\mathcal{O}\,(N^{\,3})$ (with a big constant in front of $N^{\,3}$). Thus, when $N$ is very large, the computational time may be prohibitive and the accumulation of round-off errors significant, even if an algorithm of complexity $\mathcal{O}\,(N^{\,2})$ can be used. It is well-known that for steep waves and cnoidal waves in shallow water, the number of Fourier modes needed for accurate resolutions must be very large, specially when using conformal mapping. Therefore, an algorithm with complexity lower than $\mathcal{O}\,(N^{\,2})$ is desirable to achieve these computations. It is the purpose of the present paper to describe such an algorithm with overall complexity $\mathcal{O}\,(N\log N)\,$, thus suitable for arbitrary precision computations of all waves of practical interest in arbitrary depth.

For a simple wave equation, [40] proposed an algorithm based on stabilised fixed point iterations for computing solitary waves. This algorithm is interesting because it is very simple to implement and has complexity $\mathcal{O}\,(N)\,$. However, Petviashvili’s method generally works only for equations with special nonlinear terms, typically autonomous equations with homogeneous nonlinearities [48]. Several variants have then been proposed to extend somehow the scope of these modified Petviashvili methods [1, 3, 31]. Instead of tweaking the algorithm, our approach here is to rewrite the Euler equations in a form suitable for their numerical resolution via the classical Petviashvili method. Doing so, we obtain an algorithm suitable to compute waves in arbitrary depth, in particular long waves with wavelength of tenths to millions times the mean water depth. To our knowledge, it is the first algorithm capable of computing such long waves.

The Petviashvili method works for the Euler equations if they are rewritten in the form of a Babenko equation [5]. This was successfully implemented to compute solitary gravity waves [14, 22]. Unfortunately, this algorithm does not work for periodic waves. In infinite depth, periodic waves were successfully computed with a modified Petviashvili method by [23]. Their algorithm does not work in finite depth, however. In order to overcome these drawbacks, we propose here a simple change of variable that transforms the Babenko equation into a form tractable with the classical Petviashvili method. Our algorithm is not designed for the most extreme waves: for a given wavelength, it works only for all wave-heights less than about ninety-nine percent of the maximum one. Thus, we are able to rapidly compute waves of practical interest in arbitrary depth (infinite, finite and shallow) and to arbitrary precision.

The paper is organised as follow. In Section 2, we present the physical assumptions and the mathematical definitions and notations. In Section 3, we introduce a conformal mapping and we give precise definitions of all the variables and parameters in the conformal space. In Section 4, we derive a Babenko equation written in a suitable form for a fast numerical resolution. The numerical algorithm is subsequently described in Section 5 and relevant numerical examples are provided in Section 6. Summary and perspectives are outlined in Section 8.

Definitions and notations

We consider steady two-dimensional potential flows due to surface gravity waves in constant depth $0pt\,$. The fluid is of (positive) constant density $\rho\,$, the pressure is zero at the impermeable free surface, while the seabed is fixed, horizontal and impermeable.

Let be $(x,\,y)$ a Cartesian coordinate system moving with the wave, $x$ being the horizontal coordinate and $y$ being the upward vertical one. The wave is ($2{\pi}/k$)-periodic111The fundamental wavenumber $k$ is zero for solitary and aperiodic waves. and $x\ =\ 0$ is the abscissa of a wave crest. By $y\ =-0pt\,$, $y\ =\ \eta\,(x)$ and $y\ =\ 0$ we denote, respectively, the equations of the bottom, of the free surface and of the mean water level (see Figure 1, left). The latter implies that $\left<\eta\right>\ =\ 0$ — $\left<\bullet\right>$ the Eulerian average operator over one spatial period (wavelength) — i.e.

[TABLE]

$a\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \eta\,(0)$ denotes the wave crest amplitude and $b\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ -\eta\,({\pi}/k)$ is the wave trough amplitude, so that $H\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ a\ +\ b$ is the total wave height. A wave steepness $\varepsilon$ is then classically defined as $\varepsilon\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ k\/H/2\,$. Variable $\phi$ denoting the velocity potential, the classical equations of motion are

[TABLE]

where $g\ >\ 0$ is the (constant) acceleration due to gravity and $B$ is a Bernoulli constant.

Let be $\psi\,$, $u$ and $v$ the stream function, the horizontal and vertical velocities, respectively, such that $u\ =\ \phi_{\,x}\ =\ \psi_{\,y}$ and $v\ =\ \phi_{\,y}\ =\ -\psi_{\,x}\,$. It is convenient to introduce the complex potential $f\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \phi\ +\ \mathrm{i}\,\psi$ (with $\mathrm{i}^{2}\ =\ -1$) and the complex velocity $w\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ u\ -\ \mathrm{i}\,v$ that are holomorphic functions of $z\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ x\ +\ \mathrm{i}\,y$ (i.e., $f\ =\ f\,(z)$ and $w\ =\ \mathrm{d}f/\mathrm{d}z$). The complex conjugate is denoted with a star (e.g., $z^{*}\ =\ x\ -\ \mathrm{i}\,y$), while over ‘bars’ denote the quantities written at the seabed — e.g., $\bar{z}\,(x)\ =\ x\ -\ \mathrm{i}\,0pt\,$, $\bar{\phi}\,(x)\ =\ \phi\,(x,\,y\!=\!-0pt)$ — and over ‘tildes’ denote the quantities written at the surface — e.g., $\tilde{z}\,(x)\ =\ x\ +\ \mathrm{i}\,\eta\,(x)\,$, $\tilde{\phi}\,(x)\ =\ \phi\,(x,\,y\!=\!\eta\,(x))\,$.222Note that, e.g., $\tilde{u}\ =\ \widetilde{\ \phi_{\,x}\,}\ \neq\ \tilde{\phi}_{\,x}\ =\ \tilde{u}\ +\ \eta_{\,x}\,\tilde{v}\,$. Free surface and bottom being streamlines, $\tilde{\psi}$ and $\bar{\psi}$ are constants. One can then take $\tilde{\psi}\ =\ 0$ or $\bar{\psi}\ =\ 0$ without loss of generality (gauge condition for the stream function).

The pressure field can be obtained from the Bernoulli equation

[TABLE]

where $p$ is the pressure divided by the density. At the free surface the pressure being zero (i.e., $\tilde{p}\ =\ 0$), the Bernoulli constant $B$ is defined averaging (2.2) applied at the free surface and using the condition (2.1), i.e.

[TABLE]

From the incompressibility and the irrotationality, it follows that $B$ can also be obtained from the expression at the bottom [13]

[TABLE]

and then, from the Bernoulli equation averaged at the bed, gives $\left<\bar{p}\right>\ =\ g0pt\,$. More generally, $B$ equals $u^{\,2}\ +\ v^{\,2}$ averaged along any streamline (in the frame of reference moving with the wave).

Let be $-c_{\,\textsc{s}}$ the mean flow velocity defined as

[TABLE]

Thus, $c_{\,\textsc{s}}$ is the phase velocity of the wave observed in the frame of reference without mean flow (Stokes’ second definition of phase celerity), that is also the frame where the wave impulse is zero (Appendix B). Another important quantity is the phase velocity $c_{\,\textsc{e}}$ observed in the frame of reference without mean velocity at the seabed (that is also the one where the circulation is zero, cf. Appendix B):

[TABLE]

This is Stokes’ first definition of phase celerity. Since the motion is irrotational, $c_{\,\textsc{e}}$ can be obtained averaging $u$ along any horizontal line $y\ =\ \mathrm{constant}$ (but not along the wavy surface, i.e., $c_{\,\textsc{e}}\ \neq\ -<\tilde{u}>$). Many other phase velocities can of course be defined, but $c_{\,\textsc{s}}$ and $c_{\,\textsc{e}}$ are two velocities of special interest here. Notice that $B\ =\ c_{\,\textsc{s}}^{\,2}\ =\ c_{\,\textsc{e}}^{\,2}$ in deep water and for solitary waves (see below), and that neither $c_{\,\textsc{s}}$ and $c_{\,\textsc{e}}$ are the linear phase velocity $c_{\,0}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \sqrt{(g/k)\tanh(k0pt)}$ if the wave amplitude is not zero. A discussion on the Bernoulli constant and phase velocities can be found in [11]. Further informations regarding these phase velocities and integral quantities are given here in the Appendices A, B and C.

Conformal mapping

Let be the change of independent complex variable $z\ \mapsto\ \zeta\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ (\mathrm{i}\,\tilde{\psi}\ -\ f)/c_{\,\textsc{r}}\,$, $c_{\,\textsc{r}}\ \neq\ 0$ being a velocity of reference. In practice, one could take $c_{\,\textsc{r}}\ =\ c_{\,\textsc{s}}$ or $c_{\,\textsc{r}}\ =\ c_{\,\textsc{e}}$ or $c_{\,\textsc{r}}\ =\ c_{\,0}$ or $c_{\,\textsc{r}}\ =\ B^{\,1/2}$ or $c_{\,\textsc{r}}\ =\ (c_{\,\textsc{e}}c_{\,\textsc{s}})^{\,1/2}\,$, for example, but another convenient choice can be made depending on the problem at hands. One should take $c_{\,\textsc{r}}\ >\ 0$ if the wave travels toward the increasing $x-$direction in a ‘fixed’ frame of reference, and $c_{\,\textsc{r}}\ <\ 0$ is the wave travels toward the decreasing $x-$direction. Here, without loss of generality, we consider only waves travelling toward the increasing $x-$direction.

This change of variable conformally maps the fluid fundamental domain

[TABLE]

into the rectangle (see Figure 1, right)

[TABLE]

where $\alpha\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \operatorname{Re}\,(\zeta)$ and $\beta\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \operatorname{Im}\,(\zeta)\,$. For convenience, we introduce the apparent wavenumber $\mathchar 22\relax\mkern-9.0muk$ and apparent depth $\;\mathchar 22\relax\mkern-14.0mud\/$ in the conformal plane

[TABLE]

that are generally different from the corresponding quantities in the physical plane. Note that $c_{\,\textsc{r}}$ does not appear in the expression of $\mathchar 22\relax\mkern-9.0muk\;\mathchar 22\relax\mkern-14.0mud\/\ =\ c_{\,\textsc{s}}\/c_{\,\textsc{e}}^{\,-1}\/k\/d\,$, so no choice of $c_{\,\textsc{r}}$ can enforce the equality $\mathchar 22\relax\mkern-9.0muk\;\mathchar 22\relax\mkern-14.0mud\/\ =\ kd\,$, the latter being obtained only if $c_{\,\textsc{e}}\ =\ c_{\,\textsc{s}}\,$. Conversely, the choice $c_{\,\textsc{r}}\ =\ \sqrt{c_{\,\textsc{e}}c_{\,\textsc{s}}}$ yields $\;\mathchar 22\relax\mkern-14.0mud\//\mathchar 22\relax\mkern-9.0muk\ =\ 0pt/k$ so, with this peculiar choice of $c_{\,\textsc{r}}\,$, the areas of the fundamental periods are identical in physical and conformal planes. For our numerical resolution, we found convenient to take $c_{\,\textsc{r}}\ =\ c_{\,\textsc{e}}$ (see Section 5 below).

Since $\mathrm{d}z/\mathrm{d}\zeta\ =\ z_{\,\alpha}\ =\ -\mathrm{i}\,z_{\,\beta}$ (subscripts denoting partial derivatives), we have the Cauchy–Riemann relations $x_{\,\alpha}\ =\ y_{\,\beta}$ and $x_{\,\beta}\ =\ -y_{\,\alpha}\,$, while the complex velocity and the velocity components are

[TABLE]

thence

[TABLE]

With these relations one can compute all the physical quantities of interest.

Conformal averaging operator

The Bernoulli constant, from the relation (2.3), is defined in the conformal plane by

[TABLE]

the integral being computed keeping $\beta$ constant, in particular at the free surface $\beta\ =\ 0\,$. The relation (3.2) shows that it is convenient to introduce the average operator over one period in the $\alpha-$variable ($\beta$ being kept constant)

[TABLE]

for any quantity $(\cdots)\,$. Thus, at the free surface and at the bottom we have respectively

[TABLE]

and conversely

[TABLE]

Averaged physical quantities being defined in the physical plane and not in the conformal plane, the connections between the averaging operators are useful to express these physical quantities in the conformal plane. Conversely, these connections are also useful to express averaged quantities in the conformal plane (needed for a numerical resolution) in their physical plane counterparts. For an easy reference, we give several such relations in the Appendix C.

Resolution of the conformal mapping

With the change of dependent variables

[TABLE]

the Cauchy–Riemann relations $X_{\,\alpha}\ =\ Y_{\,\beta}$ and $X_{\,\beta}\ =\ -Y_{\,\alpha}$ hold, while the bottom ($\beta\ =\ -\;\mathchar 22\relax\mkern-14.0mud\/$) and the free surface ($\beta\ =\ 0$) impermeabilities yield

[TABLE]

At the boundaries of the fundamental period (i.e., $\alpha\ =\ 0$ and $\alpha\ =\ 2\,\pi/\mathchar 22\relax\mkern-9.0muk$), we have from (3.8)(a)

[TABLE]

and more generally $X\,(\alpha\ +\ 2\pi/\mathchar 22\relax\mkern-9.0muk,\,\beta)\ =\ X\,(\alpha,\,\beta)\,$. Therefore, the function $X$ is $(2\pi/\mathchar 22\relax\mkern-9.0muk)-$periodic.

The functions $X$ and $Y$ can be expressed in term of $\bar{X}$ — i.e., the function $X$ written at the bottom — as (see [8, 9] for detailed derivations)

[TABLE]

where a star denotes the complex conjugate. Thus, the Cauchy–Riemann relations and the bottom impermeability are fulfilled identically. At the free surface $\beta\ =\ 0\,$, (3.9) yields

[TABLE]

that can be inverted as

[TABLE]

The relation (3.10) can then be rewritten with quantities expressed at the free surface only, i.e.

[TABLE]

where $\mathscr{H}$ is an anti-adjoint pseudo-differential operator acting on a pure frequency as

[TABLE]

Indeed, the differential operator $\partial_{\,\alpha}$ corresponding to a simple multiplication by $\mathrm{i}\,\kappa$ in Fourier space ($\kappa$ the frequency333$\kappa\ =\ n\,k$ for the n$${}^{\mathrm{\small\textsf{th}}} Fourier mode of a ($2\pi/k$)-periodic function.), the pseudo-differential operators are defined in Fourier space substituting $\mathrm{i}\kappa$ for $\partial_{\,\alpha}\,$. Thus, for instance, the operator $\cos\,[(\beta\ +\ \;\mathchar 22\relax\mkern-14.0mud\/)\partial_{\,\alpha}]$ is in Fourier space a multiplication by $\cosh\,[\kappa\,(\beta\ +\ \;\mathchar 22\relax\mkern-14.0mud\/)]\,$. It should be noted that the formal Taylor expansions — such as in (3.9) and (3.10) — of the pseudo-differential operators are practically useless, specially in deep water ($0pt\ \to\ \infty$) where the operator $\mathscr{H}$ becomes the classical Hilbert transform.

Averaging of dependent functions

At the bottom, averaging $x$ over one wavelength, one obtains easily

[TABLE]

From the definition of the average operator (3.3), we have also

[TABLE]

Comparing (3.13) and (3.14), we obtain

[TABLE]

meaning that $\bar{X}$ is a periodic function averaging zero. Thus, with the relations (3.11) to (3.15), the boundedness and periodicity of $\bar{X}$ imply that

[TABLE]

Hence, the functions $X$ and $Y$ have zero average in the $\alpha-$variable. The relation (3.11) can thus be inverted without ambiguities giving, in particular,

[TABLE]

$\mathscr{C}$ being a self-adjoint positive-definite pseudo-differential operator such that

[TABLE]

Note that $\mathscr{C}_{\,\infty}^{\,-1}$ is singular for $\kappa\ =\ 0\,$, but the special choice $\mathscr{C}_{\,\infty}^{\,-1}\,\left\{1\right\}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ 0$ does not matter as long as $\mathscr{C}_{\,\infty}^{\,-1}$ is applied to a function averaging to zero (or to any known value that can be enforced).

In summary, we have obtained the special relations

[TABLE]

and the averaged quantities.

[TABLE]

Mean level condition

The definition of the mean level in the transformed domain remains to be considered. Using the relations (3.8) and (3.20), the mean level condition (2.1) becomes

[TABLE]

thence

[TABLE]

This relation is fundamental, in particular for a numerical resolution. With the previous results, several useful averaged relations can be expressed in term of the important parameter \left\llangle\tilde{y}\right\rrangle (cf. Appendix C).

Celerities

The definition (3.2) of the Bernoulli constant written at the free surface and at the bottom yields

[TABLE]

the second equality deriving from the relations

[TABLE]

where $\mathscr{S}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \partial_{\,\alpha}\csc\,\left[\/\;\mathchar 22\relax\mkern-14.0mud\/\,\partial_{\,\alpha}\/\right]$ is a pseudo-differential operator acting on a pure frequency as

[TABLE]

Note that, as $|\kappa|\ \to\ \infty\,$, $\mathscr{S}$ decays exponentially fast in Fourier space, unlike $\mathscr{C}$ that grows linearly. Note also that $\mathscr{S}\,\left\{\tilde{y}\right\}\ \to\ 0$ as $0pt\ \to\ \infty\,$, hence $B\ =\ c_{\,\textsc{e}}\,c_{\,\textsc{s}}$ in deep water (together with $c_{\,\textsc{e}}\ =\ c_{\,\textsc{s}}\ =\ c$ as mentioned above) and thus $B\ =\ c^{\,2}\,$.

Babenko equations

Using (3.1)(a) and the relation $u^{\,2}\ +\ v^{\,2}\ =\ w\,w^{*}$, the Bernoulli equation (2.2) at the free surface can be written

[TABLE]

$w\ =\ u\ -\ \mathrm{i}\,v$ being a holomorphic function such that $\operatorname{Im}\,(w)\ =\ 0$ at the bottom, we have at the free surface — see the derivation of (3.11) in the previous section —

[TABLE]

thence, with (3.18) and (4.1),

[TABLE]

After simplifications and applying the antiderivative operator $\partial_{\,\alpha}^{\,-1}\,$, one obtains at once the Babenko equation

[TABLE]

where $K_{\,1}$ is an integration constant. Averaging the equation over one wavelength, then using the relations (C.3) and (C.5), one obtains an expression for the constant $K_{\,1}$

[TABLE]

For numerical resolutions, it is convenient to make the change of dependent variable $\tilde{y}\,(\alpha)\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \Upsilon\,(\alpha)\ +\ \delta\,$, where $\delta$ is a constant at our disposal. Thus, the equation (4.2) becomes

[TABLE]

where

[TABLE]

and we have (see Appendix C for details)

[TABLE]

giving two equations for $\delta$ and $c_{\,\textsc{s}}\//\/c_{\,\textsc{e}}\,$. Alternatively, a related equation is obtained obtained applying the operator $\mathscr{C}$ to (4.4)

[TABLE]

In practice, the parameter $\delta$ may be conveniently chosen such that:

• (i)

$\delta\ =\ 0\,$, so we are dealing with the free surface $\tilde{y}\,$;

• (ii)

\delta\ =\ \left\llangle\,\tilde{y}\,\right\rrangle\,, so we are dealing with the zero-mean dependant variable $\tilde{Y}\,$;

• (iii)

$\delta\ =\ \min\,(\tilde{y})\,$, so $\Upsilon\ \geqslant\ 0$ and the wave appears to somehow look like a “solitary wave” on $[\,-\pi/k,\,\pi/k\,]\,$;

• (iv)

$K_{\,2}\ =\ 0\,$, so the equation (4.4) is homogeneous.

Properties (iii) and (iv) are both desirable features for using Petviashvili’s iterations, leading to a simple and fast numerical scheme. However, no choice of $\delta$ can enforce (iii) and (iv) simultaneously, so we proceed as follow.

Note first that in deep water $0pt\ \to\ \infty\,$, $\mathscr{C}^{\,-1}\ \to\ \mathscr{C}_{\,\infty}^{\,-1}$ is singular and the equation (4.4) is problematic because $\mathscr{C}_{\,\infty}^{-1}$ is applied to the quantity $\Upsilon\,\mathscr{C}_{\,\infty}\,\left\{\Upsilon\right\}$ that has a non-zero (unknown a priori) mean value. This problem is overcome applying the operator $\mathscr{C}_{\,\infty}$ to (4.4) (the constant $K_{\,2}$ then vanishes), thus yielding the alternative homogeneous equation

[TABLE]

where

[TABLE]

The equation (4.7) for deep water with $\delta\ =\ 0$ was solved by [23] using the generalised Petviashvili method (GPM) because the classical Petviashvili method (CPM) does not converge (with $\delta\ =\ 0$). Here, we solve the equation (4.7) with $\delta\ =\ \min\,(\tilde{y})$ using classical Petviashvili’s iterations that converge in deep water and in finite depth.

Numerical resolution

We detail here how the numerical algorithm is implemented. We mostly focus on periodic waves ($k\ >\ 0$), the case of solitary waves ($k\ =\ 0$) being subsequently briefly described. Details for alternative computations of solitary waves are given in [14, 22].

User provided parameters

For practical applications, the user generally wants to define the wave by the parameters $g\,$, $0pt\,$, $k$ and $H\,$. Two of these parameters can be freely chosen to compute the solution. In deep water and finite depth (i.e. not shallow), a wave is most often defined only by the two dimensionless parameters $k0pt$ and $\varepsilon\ =\ kH/2\,$, from which the four physical parameters are defined as follow in the numerical resolution. In shallow water, the dimensionless height $H/d$ is often used instead of $\varepsilon\,$. It should be noticed that $\varepsilon$ is not a good parameter to characterise the steepness in shallow water ($\varepsilon\ =\ 0$ for all solitary waves), while $H/d$ is not suitable in deep water ($H/d\ =\ 0$ for all waves in infinite depth). Alternative definitions of the wave steepness, valid for all depths and all wavelengths, could be introduced such as $kH/\tanh(kd)\,$. However, this alternative steepness parameter being generally not used by practitioners, it is not considered here.

The water depth is considered infinite if

[TABLE]

where $\mathsf{tol}$ is the numerical tolerance used for the computations. Working with $n-$digit arithmetics, we generally take $\mathsf{tol}\ =\ 10^{\,2\,-\,n}$ yielding $\mathsf{tol}\ =\ 10^{\,-14}$ in double precision, $\mathsf{tol}\ =\ 10^{\,-32}$ in quadruple precision and $\mathsf{tol}\ =\ 10^{\,-69}$ in octuple precision (standard IEEE 754–2008). If the inequality (5.1) is fulfilled, we consider that $0pt\ =\ \infty$ and we choose (dimensionless) units such that $g\ =\ k\ =\ 1$ (this choice is always possible via a suitable scaling). Otherwise, we choose $g\ =\ 0pt\ =\ 1$ (that is also always possible without loss of generality). Many other scaling could be used, but the ones above are the most common in both deep water and finite depth.

The parameters $g\,$, $0pt$ and $k$ being now defined, and the parameter $\varepsilon$ being chosen, the total wave height is obviously $H\ =\ 2\varepsilon/k\,$. Note that the wave can be defined by other pairs of dimensionless parameters than $\{kd,\,\varepsilon\}$ or $\{kd,\,H/d\}\,$. For instance, a parameter like $|\tilde{u}_{\,0}/c_{\,\textsc{r}}|$ ($\tilde{u}_{\,0}$ the horizontal velocity at the crest) is often used to investigate the (almost) highest waves [35, 36]. However, this parameter involves quantities generally unknown from, e.g., experimental measurements, so it is not convenient for a practical use. Thus, defining a wave by such parameters would require important pre- and post-processing for practical applications, that are not always trivial.

Computational parameters

In addition to the physical parameters defined above, it is convenient (but unessential) to take $c_{\,\textsc{r}}\ =\ c_{\,\textsc{e}}$ and to introduce the dimensionless parameter $\sigma\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ c_{\,\textsc{s}}/c_{\,\textsc{e}}$ (hence $\mathchar 22\relax\mkern-9.0muk\ =\ k\,$, $\;\mathchar 22\relax\mkern-14.0mud\/\ =\ \sigma 0pt$). In infinite depth $\sigma\ =\ 1$, but in finite depth $\sigma$ is unknown and therefore must be computed (see below). We also chose $\delta\ =\ \min(\tilde{y})$ because we shall use Petviashvili’s iterations.

The fundamental period is discretised with $2\,N$ equally spaced points $\mathchar 22\relax\mkern-9.0muk\/\alpha_{\,j}\ =\ (j\ -\ 1)\pi/N$ for $j\ =\ 1,\,2,\,\ldots,\,2\,N\,$. The number of nodes should be large enough to ensure the desired accuracy. For infinitesimal short waves (i.e. $\varepsilon\ \ll\ 1$ and $k\,d$ not small), $N\ =\ 32$ (say) may be sufficient to ensure machine precision when working in double precision. However, the number of required nodes increases rapidly as the wave-height and the wavelength increase. Thus, for steep and cnoidal waves, the number of nodes must be substantially larger. Since we are solving the equation with an algorithm requiring $\mathcal{O}\,(N\log N)$ operations, $N$ can be large without significant speed burden. Therefore, one can take $N\ =\ 1\,024$ as default when working in double precision because it is often sufficient. If more nodes are needed, the cost of increasing $N$ is low because the algorithm is very fast. However, $N$ cannot be too large in order to avoid large accumulation of round-off errors (when working in double precision, it is wise to take $N\ \leqslant\ 2^{\,19}\ =\ 524\,288\,$, say). If an extremely large number of nodes is required, then quadruple or higher precision is necessary to avoid a significant accumulation of round-off errors.

Classical Petviashvili method (CPM)

Petviashvili’s iterations are a simple and efficient algorithm for computing solitary waves solution of the Babenko equation [14, 22]. This algorithm is very easy to implement and runs fast since each iteration requires $\mathcal{O}\,(N)$ operations. Note that Petviashvili’s iterations are used here together with fast Fourier transforms (FFT) that require $\mathcal{O}\,(N\log N)$ operations, so the overall complexity of our algorithm is $N\log N\,$.

In order to apply Petviashvili’s method, we separate the linear and nonlinear terms and rewrite the equation (4.7) as $\mathscr{L}\,\{\Upsilon\}\ =\ \mathscr{N}\,\{\Upsilon\}$ with

[TABLE]

where $\sigma\ =\ c_{\,\textsc{s}}/c_{\,\textsc{e}}\,$, $\delta\ =\ \min(\tilde{y})$ and $\Upsilon\ =\ y\ -\ \delta\,$. From the approximation $\Upsilon^{\,(i)}$ obtained at the i$${}^{\mathrm{\small\textsf{th}}} iterations, the next approximation $\Upsilon^{\,(i+1)}$ is obtained via the Petviashvili weighted fixed point iteration

[TABLE]

At each iteration, the zero value of $\Upsilon^{\,(i)}$ at the trough and the wave height are enforced via the renormalisation

[TABLE]

This renormalisation improves the convergence, specially for steep waves.

In finite depth, the operator $\mathscr{L}^{\,-1}$ is singular but $\mathscr{L}^{\,-1}\circ\mathscr{N}$ is regular. This is because both operators involve the factor $\mathscr{C}_{\infty}$ that is zero for the zero frequency (this factor is introduced to kill the constant $K_{\,2}$). Doing so, the Babenko equation has been singularised at the zero frequency, but this singularity is only apparent (i.e. movable in the sense that $\mathscr{C}_{\infty}^{\,-1}\circ\mathscr{C}_{\infty}$ is identity). However, defining explicitly $\mathscr{L}^{\,-1}\circ\mathscr{N}$ at the zero frequency is not necessary, the mean value of $\Upsilon^{\,(i+1)}$ being enforced by the renormalisation (5.5). Indeed, setting arbitrarily $\mathscr{L}^{\,-1}\,\{1\}\ \equiv\ 0\,$, $\Upsilon^{\,(i+1)}$ computed with (5.4) is obtained modulo an unknown (generally incorrect) mean value, the right value of \left\llangle\Upsilon\right\rrangle being subsequently enforced via (5.5).

As initial guess $\Upsilon^{\,(0)}\,$, we take the linear approximation

[TABLE]

unless a better guess is provided by the user, for instance from another calculation with slightly different parameters $k0pt$ and $\varepsilon$ (useful for analytic continuations). We found that, with the initial guess (5.6), the Petviashvili iterations always converge, even for large waves in shallow water, so we did not try out other guesses. The convergence from the initial guess (5.6) illustrates the robustness of the method.

Petviasvili’s iterations (5.4) involve the unknown parameters $B\,$, $\sigma$ and $\delta$ via the definition of the operators $\mathscr{C}\,$, $\mathscr{L}$ and $\mathscr{N}\,$. Therefore, $B\,$, $\sigma$ and $\delta$ must be computed from $\Upsilon^{\,(i)}$ before (5.4) can be used. (It would be the same with any other iterations, such as Newton and Levenberg–Marquardt methods.) These parameters are obtained as follow.

Computation of the unknown parameters

We first compute \tilde{Y}\ =\ \Upsilon^{\,(i)}\ -\ \left\llangle\,\Upsilon^{\,(i)}\,\right\rrangle\,. In deep water $\sigma\ =\ 1$ and, according to (3.21), we have

[TABLE]

thence \tilde{y}\ =\ \tilde{Y}\ +\ \left\llangle\,\tilde{y}\,\right\rrangle by definition of $\tilde{Y}\,$.

In finite depth, in general $\sigma\ \neq\ 1$ is unknown ($\sigma\ =\ 1$ only for solitary waves) and must be computed. To do so, the relation (3.21) is rewritten as the equation

[TABLE]

It should be recalled here that $\mathscr{C}$ depends on $\sigma\,$, so (5.8) is a nonlinear equation for $\sigma\,$. Equation (5.8) is thus solved with Newton iterations

[TABLE]

with $\mathscr{S}\ =\ \partial_{\alpha}\csc\!\left[\/\sigma\/0pt\/\partial_{\alpha}\/\right]\,$. In practice, one Newton iteration is sufficient because the initial guess $\sigma_{\,0}$ is given by the approximation of $\sigma$ obtained at the previous iteration from $\Upsilon^{\,(i-1)}$ that, if the CPM converges, is closed to the exact solution. Once $\sigma$ has been obtained, we compute \left\llangle\tilde{y}\right\rrangle\ =\ (\sigma\,-\,1)\,d and \delta\ =\ \left\llangle\tilde{y}\right\rrangle\ -\ \left\llangle\Upsilon\right\rrangle\,, so these parameters are now known for the i$${}^{\mathrm{\small\textsf{th}}} Petviashvili iteration. It should be emphasised that accurate computations of $\sigma$ and \left\llangle\tilde{y}\right\rrangle are absolutely crucial to ensure that the mapping $z\ \mapsto\ \zeta$ is conformal and that the still water level is where it should be.

Finally, the Bernoulli constant $B$ is obtained from the equation (4.6) applied at the crest ($\alpha\ =\ 0$) and at the trough ($\alpha\ =\ \pi/\mathchar 22\relax\mkern-9.0muk$), i.e.,

[TABLE]

thence

[TABLE]

with the notations $[f]_{\,a}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ f(a)$ and $[f]_{\,a}^{\,b}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ f\,(b)\ -\ f\,(a)\,$. Relation (5.12) is obtained subtracting (5.10) and (5.11), thus $K_{\,2}$ vanishes and does not need to be computed.

All the parameters involved in the Babenko equation are now defined and the Petviashvili iterations (5.4) can be applied until the desired accuracy is reached.

Post processing

After convergence of the Petviashvili iterations, all the parameters of interest can be computed. The celerity $c_{\,\textsc{e}}$ is given by

[TABLE]

thence $c_{\,\textsc{s}}\ =\ \sigma c_{\,\textsc{e}}\,$. These parameters being defined, all the integral quantities in Appendix B are easily computed.

Often, users want to know the velocity and other fields inside the bulk at a given location $z\,$. This can be obtained from the integrals provided in the Appendix D and discretised according to the trapezoidal rule [22]. This is very simple to implement and also very accurate, provided that $z$ is not too close to the free surface. (Typically, the distance between $z$ and the free surface should be larger than $\Delta\alpha\,$.)

Remark 1**.**

Since we use a spectral method to solve a nonlinear pseudo-differential equation, some aliasing errors may occur. With a number of Fourier modes large enough to achieve spectral accuracy, the magnitude of all the neglected high-frequency modes is smaller than the machine precision (about $10^{\,-16}$ in double precision). The equation being quadratic in nonlinearity, in theory a two-thirds rule [6] should be applied in order to avoid aliasing. However, even if the two-thirds rule is not applied, the aliasing error is about machine precision, provided that the aliased frequencies (i.e. the upper third of the spectra) are small (say of magnitudes of at most $10^{\,-8}$ in double precision) so that their products are numerically zero. This means that an anti-aliasing filter is not necessary if $N$ is large enough, as can be seen in the spectra below. Of course, one can easily introduce an anti-aliasing (i.e. low-pass) filter if this turns out to be necessary for some computations.

Numerical examples

The algorithm described above has been implemented in MatlabTM and it is freely downloadable [15]. This program was written with clarity in mind, so it can be easily understood, modified and translated into any programming language. In particular, this program can be easily modified to run in arbitrary precision, provided that this feature is available to the user (we use the Advanpix Multiprecision Toolbox [27]).

Deep water

In deep water ($d\ \to\ \infty$), periodic waves with identical crests are obtained with our algorithm provided that $\varepsilon\ \lesssim\ 0.44\,$. Thus, the algorithm converges for rather steep waves (up to about $99.3\%$ of the highest waves), the maximum steepness being $\varepsilon\ \approx\ 0.443164$ [35].

It should be noted that with $\delta\ =\ 0$ the CPM diverges for all steepnesses. This shows that the choice $\delta\ =\ \min(\tilde{y})$ improves significantly the convergence of the CPM, rending unnecessary the use of the GPM. It should also be noted that Fenton’s algorithm [24] converges for $\varepsilon\ \lesssim\ 0.36$ and that it is much slower than the method described here. For $\varepsilon\ \lesssim\ 0.36\,$, Fenton’s and ours algorithms match up to about six digits (for non-infinitesimal waves), that corresponds to the accuracy of Fenton’s algorithm.

Actually, the present algorithm can converge also for $\varepsilon\ >\ 0.44\,$, but then ‘ghost’ solutions [21] are obtained (similar to the one on the right Figure 2). For $\varepsilon\ =\ 0.44\,$, varying $N$ leads to the same solution in double and quadruple precisions. However, for $\varepsilon\ =\ 0.4401\,$, with $N\ =\ 1\,024$ a ‘normal’ solution is obtained (Figure 2 left), but a ‘ghost’ (spurious) solution is obtained with $N\ =\ 2\,048$ (Figure 2 right). The algorithm behaviour for the highest computable waves is discussed in the Section 7 below.

For $\varepsilon\ \leqslant\ 0.44$ the algorithm converges rapidly444On a 2012 MacBook Pro laptop computer, with $N\ =\ 2\,048\,$, $\mathsf{tol}\ =\ 10^{\,-15}\,$, $\varepsilon\ =\ 0.4$ and in double precision, the solution is obtained in about half a second. to the solution. Actually, any arbitrary accuracy can be achieved provided that $N$ is large enough (Figure 3). For instance, for $\varepsilon\ =\ 0.4\,$, the solution is obtained to machine double precision with $N\ =\ 512$ (Fig. 3 upper). $N\ =\ 1\,024$ is not sufficient to achieve machine quadruple precision (Fig. 3 middle), the latter being obtained for $N\ =\ 2\,048$ (Fig. 3 lower). However, $N\ =\ 2\,048$ is not sufficient to achieve full octuple precision, that can be obtained with larger $N\,$. Similarly, any accuracy can be obtained provided that $N$ is large enough. Note that the Figure 3 clearly shows that aliasing errors are not significant although no special anti-aliasing techniques were applied.

This test illustrates the accuracy and the robustness of the algorithm. Indeed, some algorithms diverge when $N$ is too large although they converge for smaller $N$ (for the same steepness); a problem not faced by the algorithm described here.

Finite depth

We found that the CPM and GPM are both divergent when applied to the Babenko equation (4.7) with $\delta\ =\ 0\,$. Conversely, the CPM converges well if one takes $\delta\ =\ \min(\tilde{y})\,$, so we did not try the GPM. As for the deep water case, the algorithm converges for all but the highest waves, i.e., the algorithm converges up to the first maximum of $B\,(\varepsilon)\,$. For higher waves, the algorithm diverges or converges to a ghost solution, the maximum steepness computable depending on the depth (see section 7).

For example, consider the case $kd\ =\ 1$ with steepnesses $\varepsilon\ =\ 0.3$ and $\varepsilon\ =\ 0.31\,$. Although these two large steepnesses are close, they correspond to quite different free surfaces (Fig. 4 upper), as can be seen in their spectra (Fig. 4 lower). With $N\ =\ 2\,048$, the case $\varepsilon\ =\ 0.3$ is resolved to machine double precision (Fig. 4 lower left), while the case $\varepsilon\ =\ 0.31$ is resolved only to a mild accuracy (Fig. 4 lower right). Machine double precision is achieved for $\varepsilon\ =\ 0.31$ with $N\ =\ 8\,192\,$, however. This shows that as the wave approches the highest one, the number $N$ of Fourier modes has to be increased dramatically (exponentially fast) in order to reach a full spectral precision. This is inherent to the solution formulated in the conformal plane and it has nothing to do with the numerical algorithm described in the present paper. Indeed, for the highest wave with an angular crest, the Fourier spectrum decays algebraically while smaller waves have spectra decaying exponentially fast at high frequencies.

This example illustrates, like the one in deep water, the need for a large number of Fourier modes in order to achieve the full precision of a given floating point format. This example also illustrates the rapid increases of $N$ necessary for an accurate resolution as the steepness increases and, therefore, the need for a fast algorithm. Since the algorithm described here has an overall complexity $\mathcal{O}\,(N\log N)\,$, the necessity of large $N$ is not problematic and accurate computations are rapidly achieved.

For $kd\ =\ 1$ the highest computable wave has steepness $\varepsilon\ \approx\ 0.3146\,$, that is approximately $99.6\%$ of the highest wave. Limiting values for other $kd$ are given in the Table 3 showing that the highest computable wave, for any given depth and wavelength, is about $99\%$ of the maximum one.

Shallow water

When the depth $d$ over wavelength $L\ =\ 2\pi/k$ ratio is very small, i.e. $kd\ll\ 1\,$, we are dealing with the so-called shallow water situation. It is well-known that Stokes’ expansion fails to approximate such solutions of finite amplitude. Some shallow water approximations have then be proposed to approximate periodic waves, the so-called cnoidal waves. Many numerical algorithms devoted to the resolution of the full equations also fail in shallow water. For instance, Fenton’s algorithm [24] works only for $L/d\ \lesssim\ 30\,$.555For $L/d\ >\ 30\,$, Fenton’s algorithm does not converge or converge to ghost solutions with spurious oscillations, see Fig. 3-2 in [26].

The present algorithm works in shallow water without difficulties. For instance, with $L/d\ =\ 1000$ and $H/d\ =\ 0.4$ (i.e., $\varepsilon\ \approx\ 0.00126$ and $kd\ \approx\ 0.00628$) the solution is obtained (in about $0.4\,\mathsf{s}$ on a MacBook Pro laptop computer from 2012) using $N\ =\ 8\,192$ Fourier modes that are necessary to achieve machine double precision (Figure 5). Another more extreme example is the case $L/d\ =\ 10000$ and $H/d\ =\ 0.7\,$, that is computed in double precision with $N\ =\ 2^{\,19}$ and $\mathsf{tol}\ =\ 10^{\,-14}$ after $237$ iterations in about $1\,\mathsf{mn}$. It should be noted that these solutions are obtained from the initial guess (5.6) that is not at all a decent approximation of the solution, thus illustrating the robustness of the algorithm.

As the steepness or the wavelength increases, the number of Fourier modes $N$ required to achieve spectral accuracy increases rapidly. For example, with $L/d\ =\ 71$ and $H/d\ =\ 0.802$ (i.e., $\varepsilon\ \approx\ 0.0355$ and $kd\ \approx\ 0.0885$) a steep cnoidal wave is obtained (Figure 6 upper), its computation to full spectral accuracy requiring $N\ =\ 2^{\,17}\ =\ 131\,072$ Fourier modes (Figure 6 middle).666This result was obtained in less than three minutes. A rough estimate suggests that the same computation with Newton or Levenberg–Marquardt iterations, instead of Petviashvili’s ones, would take several days (possibly weeks) on the same computer.

For steep cnoidal waves in very shallow water the number of necessary Fourier modes can be prohibitively large. An alternative is to compute a shorter cnoidal wave and to eventually increase the length of the trough. Indeed, a cnoidal wave surface decaying rapidly from the crest, it rapidly reaches its minimum to machine precision (Figure 6 lower). The example of Figure 6 shows that longer cnoidal waves can be obtained to machine precision increasing the length of the trough up to the desired wavelength, then redefining the mean water level and the mean depth, as well as all the related parameters (renormalisation). A similar procedure can be used to compute solitary waves, as shown below.

Comparison with KdV cnoidal wave

Korteweg and de Vries [30] proposed an analytic approximation for small amplitude long periodic waves in shallow water. They coin the term ‘cnoidal’ wave because this approximation can be expressed in term of the Jacobi cn$-$function. In our notations, KdV analytic solution can be conveniently written

[TABLE]

dn being the elliptic functions of Jacobi of parameter $m$ ($0\ \leqslant\ m\ \leqslant\ 1$), $K\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \text{K}\/(m)$ and $E\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \text{E}\/(m)$ being the complete elliptic integrals of the first and second kinds, respectively [2].

Though an analytic approximation, KdV cnoidal wave requires significant computations. Indeed, a wave being generally defined for given height $H$ and wavenumber $k$, the parameter $m$ must be determined solving numerically the equations in (6.1) relating the parameters. For very long waves, $m$ is very close to one, to an extend that it cannot be practically computed (see Tables 1 and 2). For instance, for the very long small amplitude cnoidal wave with $L/d\ =\ 10^{\,6}$ and $H/d\ =\ 10^{\,-1}$ we have $1\ -\ m\ \approx\ 5.78\times 10^{\,-118936}\,$, a value that cannot be easily computed. This problem becomes more severe as the amplitude increases (Table 2). Thus, for very long waves, KdV analytic cnoidal solution is useless for practical applications, even if only crude approximations are sufficient. It is then more efficient to solve numerically the KdV equation, for instance with FFT and Petviashvili’s method as illustrated here for the Babenko equation. But, doing so, solving KdV is not much less demanding than solving Babenko, so the latter should be preferred. For the extreme example777This example is given to illustrate the efficiency of the method, not for its practical interest. $H/d\ =\ 10^{\,-1}$ and $L/d\ =\ 10^{\,6}\,$, Babenko solution is computed to double-precision spectral precision with $N\ =\ 2^{\,21}$ in about $100\,\mathsf{s}$ using our algorithm,888On the same hardware, for the same wave, rough estimates indicate that algorithms of complexity $\mathcal{O}\,(N^{\,2})$ would take months, and those of complexity $\mathcal{O}\,(N^{\,3})$ would take thousands of years. while KdV analytic solution cannot be computed in double precision.999For a direct determination of $m$ from the relations (6.1), one would have to use something like two hundred thousand digits computation. An alternative is to derive better conditioned relations via some non trivial mathematical manipulations, thus loosing the analytical simplicity of KdV. Of course, this drawback is not limited to the KdV analytic solution, it is a fortiori present in all cnoidal-like approximations, such as the solutions of the Boussinesq-like equations. These considerations demonstrates that our algorithm for solving the irrotational Euler equations is also a suitable alternative to simple analytic models.

Solitary waves

Solitary waves decaying exponentially fast, their surface elevation reaches zero to machine precision close to the crest. Thus, solitary waves can be efficiently computed in a periodic box, provided that the box is long enough so the periodisation does not affect the solution. This numerical trick is well-known and has been used by many authors.

The steep cnoidal example of Figure 6 reaches its minimal elevation to machine precision before its trough at $x\ =\ L/2$ (see Figure 6 lower). Thus, this cnoidal wave can be considered as a solitary wave computed in a periodic box, but with a different still water level. The actual depth for the solitary wave is $d_{\infty}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ d\ +\ \eta\,(L/2)\ =\ d\ -\ b\,$, the surface elevation is $\eta_{\infty}\ =\ \eta\ +\ b$ and the dimensionless amplitude is $H/d_{\infty}\,$. For the example of Figure 6 we obtain $H/d_{\infty}\ =\ 0.8236847804878956\,$. This result is surprising because the direct computation of solitary waves via the CPM converges only for $H/d\ \lesssim\ 0.79$ [22, 14]. The solitary wave thus obtained has been compared to an approximation obtained with Tanaka’s method using $1\,024$ nodes (Tanaka’s algorithm [45] is way too slow to use it with $2^{\,18}$ nodes). We found that the two solutions match up to about five digits, that is consistent with the accuracy of Tanaka’s method [22, 14], confirming the solution obtained by the CPM after renormalisation.

Remarks on the highest computable waves

In order to characterise the highest computable waves, we consider several dimensionless parameters. For a given $kd\,$, we denote $\hat{\varepsilon}\,$, $\varepsilon_{\mbox{\tiny\sc b}}$ and $\varepsilon_{\mbox{\tiny\sc c}}$ the steepnesses corresponding of, respectively, the highest wave (with a $120^{\circ}$ inner angle at the crest), the first maximum of the Bernoulli constant $B$ and the highest computable wave. We also consider the maximum inclination of the free surface $\theta_{\,\text{max}}$ and the parameters

[TABLE]

where $\tilde{u}_{\,0}$ and $\tilde{u}_{\,1}$ are the horizontal velocities at the crest and at the trough, respectively. Note that $\vartheta\ \to\ \infty$ and $A\ \to\ \infty$ as $\varepsilon\ \to\ \hat{\varepsilon}\,$, a wave is then considered steep if these parameters are larger than two [35].

The maximum computable steepness $\varepsilon_{\mbox{\tiny\sc c}}$ is determined (by dichotomy) for each $kd$ only up to the fourth decimal place, that is sufficient here for the discussion. All the computations of $\varepsilon_{\mbox{\tiny\sc c}}$ in Table 3 were performed in double precision, with $N\ =\ 131\,072$ Fourier positive modes and with tolerance $\mathsf{tol}\ =\ 10^{\,-12}$ for the iterations. The corresponding parameters $A\,$, $\vartheta$ and $\theta_{\text{max}}$ are given by truncated decimal expansions (i.e., not rounded to the nearest decimal approximation). In the Table 3, the steepnesses of the highest waves $\hat{\varepsilon}$ and the steepnesses $\varepsilon_{\mbox{\tiny\sc b}}$ corresponding to the first maximum of the Bernoulli constants were kindly provided by Professor Dmitri Maklakov who guaranteed that the six decimals are correct. The values of $\varepsilon_{\mbox{\tiny\sc b}}$ were confirmed by our algorithm, thus providing another validation of the method.

As shown in Table 3, the highest computable waves are rather steep with $\theta_{\text{max}}\ \approx\ 30^{\circ}\,$, $A$ and $\vartheta$ being significantly larger than two. More interestingly, $\varepsilon_{\mbox{\tiny\sc c}}$ always exceeds $\varepsilon_{\mbox{\tiny\sc b}}\,$, $\varepsilon_{\mbox{\tiny\sc c}}$ corresponding at least to $99\%$ of the maximum steepness $\hat{\varepsilon}\,$. This is somewhat surprising because similar algorithms for periodic waves in deep water [23] and solitary waves [22, 14] have maximum computable steepnesses $\varepsilon_{\mbox{\tiny\sc c}}\ \approx\ \varepsilon_{\mbox{\tiny\sc b}}\,$.

Indeed, in their implementation of the GPM for the Babenko equation (4.7) with $\delta\ =\ 0\,$, [23] found convergence for $\varepsilon\lesssim 0.436$, and for solitary waves the present authors [22, 14] noticed that the CPM converges for $H/d\ \lesssim\ 0.79\,$. These maximum computable steepnesses correspond to the first maximum of the Bernoulli constant (Figure 7, Table 3). The limiting value $\varepsilon_{\mbox{\tiny\sc c}}\ \approx\ \varepsilon_{\mbox{\tiny\sc b}}$ was then explained by the fact that, for steeper waves, the lack of one-to-one correspondance between the parameters prevents the algorithm to ‘decide’ which solution should be retained.

The situation is a bit different here, where we found that systematically $\varepsilon_{\mbox{\tiny\sc c}}\ >\ \varepsilon_{\mbox{\tiny\sc b}}$ (see Table 3), even for solitary waves (see Section 6.5). Clearly, the modifications and rescaling, introduced to make the CPM work for arbitrary wavelengths, have a beneficial effect on the steepest computable waves. A possible explanation is that the rescaling somehow enlarge the region of one-to-one correspondance between the computed parameters. This is a conjecture requiring rigorous mathematical investigations that are far beyond the purpose of the present paper.

Discussion

We described an efficient algorithm for computing steady surface gravity waves for an ideal homogeneous fluid in irrotational motion. After analytic transformations (conformal mapping, rewriting of the conditions at the free surface, change of dependent variables) we ended up with a modified Babenko-like equation that can be solved numerically via the classical Petviashvili method (CPM). The algorithm thus obtained is very fast — with complexity $\mathcal{O}\,(N\log N)\,$, $N$ being the number of Fourier modes — and any accuracy can be reached for all depths, provided that $N$ is large enough and that the steepness is not too close to the limiting one. All waves of practical interest can therefore be computed. To our knowledge, it is the first algorithm uniformly valid for all wavelength-over-depth ratios, in practice and not only in principle, that is moreover accurate in the sense that arbitrary precision numerical solutions can be obtained.

With this algorithm, the computation of steady gravity waves for the irrotational Euler equation is not more demanding than the numerical resolution of simplified models such as KdV. It should be noted that some simplified water wave models — e.g., some variants of the Boussinesq equations for long waves in shallow water — have inhomogeneous nonlinear term. Therefore, the (like of) Petviashvili method does not work and Newton or Levenberg–Marquardt methods should be used instead. The latter having complexity $\mathcal{O}\,(N^{\,3})\,$, the numerical resolution of these simplified models is more demanding than our algorithm for the irrotational Euler equations.

It is often believed that the CPM works only for localised (solitary) waves. Here, we disproved this belief with strong numerical evidences. Rigorous mathematical results are scarce [4, 39] for the (like of) Petviashvili method. Deeper mathematical understanding would be beneficial for improving the method, in particular for the computation of almost highest waves. We hope that the numerical evidences presented here will stimulate such investigations.

There exist steady surface gravity waves with different crests and asymmetric waves [38]. However, these solutions exist close to the limiting ones where the classical Petviashvili method does not work, at least not as formulated in this paper. Improving this fixed point iteration method may increase the range of computable steepnesses. However, this will not be sufficient for the computation of extreme waves using Fourier decomposition together with the conformal mapping, because a huge number of Fourier modes is needed to achieve high accuracy. This can be understood considering the limiting wave with an $120^{\circ}$ inner angular crest. Such solutions have a power $2/3$ singularity at the crest in the conformal plane [43] and, therefore, the $n$-th Fourier coefficient decays like $n^{-5/3}$ as $n\ \to\ \infty\,$. Thus, truncating the Fourier expansion after the N$${}^{\mathrm{\small\textsf{th}}} term, the error decays like $N^{\,-2/3}$ as $N\ \to\ \infty\,$. In practice, it means that in order to increase the accuracy by two digits, the number of computed Fourier modes must be multiplied by (roughly) one thousand! Clearly, the ‘brutal force’ approach consisting in massively increasing $N$ is inefficient for extreme waves, even on a powerful computer. To overcome this difficulty, several authors have proposed a change of independent variable such that the corresponding Fourier spectrum decays faster. Another possibility for accurate computations of extreme waves consists in a totally different mathematical reformulation of the problem, so that the numerical resolution is less demanding, as suggested in [12]. Needless to say, such approaches are suitable if the number of modes $N$ required for high accuracy is relatively small if the algorithmic complexity is $\mathcal{O}\,(N^{\,2})$ or $\mathcal{O}\,(N^{\,3})\,$, or if the algorithmic complexity remains $\mathcal{O}\,(N\log N)$ with $N$ not too large. To our knowledge, an efficient algorithm for arbitrary precision calculation of extreme waves has yet to be discovered.

Here, we introduced simple tricks in order to successfully apply the classical Petviashvili method for computing steady water waves of arbitrary wavelength, in arbitrary depth and to arbitrary precision. These tricks can certainly be used for other equations involved in fluid mechanics.

Acknowledgments

This work was supported by the Spanish Ministerio de Economía y Competitividad under the Research Grant MTM2014-54710-P. The authors would like to thank Professor Angel Durán for useful discussions. The authors are also grateful to Professor Dmitri Maklakov for providing some valuable data in the Table 3 and whose comments helped to improve the paper.

Appendix A Potential and stream function in a ‘fix’ frame of reference

The definition (2.6) of $c_{\,\textsc{e}}$ implies that the function $\Phi\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \phi\ +\ c_{\,\textsc{e}}\/x$ averages zero along any horizontal line $y\ =\ \text{constant}\,$, in particular at the bed $y\ =\ -0pt\,$. This interesting property suggests the introduction of a stream function $\Psi$ in this frame of reference where the mean horizontal velocity is zero at the bed and such that

[TABLE]

thence

[TABLE]

The definitions of $c_{\,\textsc{e}}$ and $c_{\,\textsc{s}}$ also yield (integrating by parts and exploiting the irrotationality)

[TABLE]

thence

[TABLE]

that is exploited in the relation (C.1). It should be noticed that the quantity $(c_{\,\textsc{s}}\ -\ c_{\,\textsc{e}})\/0pt$ (related to the wave impulse, see Appendix B) must remain bounded as $0pt\to\infty$, implying that $c_{\,\textsc{s}}\to c_{\,\textsc{e}}$ as $0pt\to\infty\,$. Note also that, the velocity magnitude varying monotonically from the bottom to the surface, we necessarily have $c_{\,\textsc{e}}\ \geqslant\ c_{\,\textsc{s}}$ [19] and therefore \left\llangle\tilde{y}\right\rrangle\ \leqslant\ 0 (see Eq. C.1).

Appendix B Integral quantities

The wave can be characterised by several integral parameters [32, 33, 37, 42]. In the frame of reference moving with the wave, there are three physically important constants: the fluid flow $\mathcal{Q}_{\,0}\,$, the momentum flux $\mathcal{S}_{\,0}$ and the energy flux $\mathcal{F}_{\,0}\,$, defined by

[TABLE]

These quantities are related to other averaged quantities of physical interest (see below). In particular, averaging $\mathcal{S}_{\,0}$ over one wavelength and exploiting the relation (2.4) and the impermeability of the free surface, one obtains at once

[TABLE]

where $\mathcal{V}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \left<\,{\textstyle{1\over 2}}\,g\,\eta^{\,2}\,\right>$ is the potential energy of the gravity force.

Other integral quantities can be defined relatively to the uniform flow of speed $-c_{\,\textsc{r}}\,$, i.e., in the fixed frame of reference where the phase velocity is $c_{\,\textsc{r}}\,$. The integral quantities of interest here are the:

[TABLE]

The equalities in these integral relations are easily obtained via some trivial derivations. Note that the radiation stress defined here differs from the definition of [32], that is $\mathcal{S}_{\,x\,x}^{\,\text{\tiny\sc LH}}\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\mathcal{S}_{\,x\,x}\ -\ \mathcal{V}\,$. Note also that the group celerity defined above is not the linear one, i.e. $c_{g}\ \neq\ \partial(k\,c_{\,0})/\partial k$ if $H\ \neq\ 0\,$.

Appendix C Relations between averaged quantities in physical and conformal planes

Averaging in the physical plane is different than averaging in the conformal plane. For practical applications, many averaged quantities may have to be computed. For an easy reference, we give below various connections between averaged quantities in the physical and conformal planes. In particular, we have the special “conformal averaged” relations at the free surface

[TABLE]

with $\tilde{y}\,(\alpha)\ =\ \eta\,(\tilde{x}\,(\alpha))$ and at the bottom

[TABLE]

with $\bar{y}=-0pt$. We also have the special “physical averaged” relations

[TABLE]

All these relations can be easily obtained from their definition and using the transformations (3.4)–(3.7) between the averaging operators. Other relations can be similarly obtained.

Appendix D Velocity and pressure fields in the fluid

In the numerical procedure described below, we use conformal mapping and a Fourier pseudo-spectral method to solve the equations. This means that we obtain a discrete approximation equally spaced along each streamline. However, for practical applications, it is often necessary to determine the fields (velocity, pressure, etc.) at various positions that are not necessarily the nodes used for the computation. These informations can be obtained as follows.

Let be $W\,(z)\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ c_{\,\textsc{e}}\ +\ w\,(z)$ the complex velocity observed in the frame of reference where the fluid velocity averages to zero at the bottom. In this frame of reference, it follows that the complex potential $F\,(z)\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}c_{\,\textsc{e}}\,z\ +\ f(z)$ (i.e., $W\ =\ \mathrm{d}F/\mathrm{d}z$) is a periodic function, bounded in the whole fluid domain ($F$ is unbounded in any other frame of reference).

The complex velocity being known at the free surface from our approximation procedure, $W$ at any complex abscissa $z$ can be obtained from the Cauchy integral

[TABLE]

where $\theta\ =\ 2\pi$ if $z$ is strictly inside the fluid domain (i.e., $\operatorname{Im}(z)\ <\ \eta$), $\theta\ =\ \pi$ if $z$ is at the free surface (i.e., $\operatorname{Im}(z)\ =\ \eta$) and $\theta\ =\ 0$ if $z$ is strictly above the free surface (i.e., $\operatorname{Im}(z)\ >\ \eta$). The bottom impermeability being taken into account via the method of images, the Cauchy integral (D.1) yields for any $z$ below the free surface

[TABLE]

where $\tilde{z}\,(\alpha)\ =\ (c_{\,\textsc{r}}/c_{\,\textsc{e}})\alpha\ +\ \tilde{X}(\alpha)\ +\ \mathrm{i}\eta(\alpha)$ and $\tilde{z}^{\,\prime}(\alpha)\ =\ \mathrm{d}\tilde{z}/\mathrm{d}\alpha\ =\ (c_{\,\textsc{r}}/c_{\,\textsc{s}})\ +\ \mathscr{C}\{\eta\}\,(\alpha)\ +\ \mathrm{i}\eta_{\alpha}(\alpha)\,$, $\tilde{X}$ and $\eta$ being known from the numerical resolution of the Babenko equation.

The integral relation (D.2) is not suitable for periodic domains. For the latter, the infinite integral is replaced by one over one period involving the Hilbert kernel

[TABLE]

From this relation, we obtain the derivative of $W$ (required to compute the acceleration field)

[TABLE]

and the complex potential

[TABLE]

such that $W\ =\ \mathrm{d}F/\mathrm{d}z$ and $\operatorname{Im}(F)\ =\ 0$ at the bed.