Stokes waves with constant vorticity: II. folds, gaps and fluid bubbles

TL;DR

This paper investigates complex wave behaviors in flows with constant vorticity, revealing folds, gaps, and fluid bubbles, and introduces efficient numerical methods to analyze these phenomena.

Contribution

It formulates the Stokes wave problem with vorticity as a nonlinear pseudodifferential equation and explores the intricate wave structures and limiting behaviors.

Findings

Folds and gaps appear in wave speed-amplitude relations for strong vorticity.

Touching waves tend to limiting Crapper waves or circular vortex waves as vorticity increases.

Multiple circular bubbles of fluid are observed in higher-gap wave profiles.

Abstract

The Stokes wave problem in a constant vorticity flow is formulated, by virtue of conformal mapping techniques, as a nonlinear pseudodifferential equation, involving the periodic Hilbert transform, which becomes the Babenko equation in the irrotational flow setting. The associated linearized operator is self-adjoint, whereby the modified Babenko equation is efficiently solved by means of the Newton-Conjugate Gradient method. For strong positive vorticity, a `fold' appears in the wave speed versus amplitude plane, and a `gap' as the vorticity strength increases, bounded by two touching waves, whose profile contacts with itself at the trough line, enclosing a bubble of air. More folds and gaps follow for stronger vorticity. Touching waves at the beginnings of the lowest gaps tend to the limiting Crapper wave as the vorticity strength increases indefinitely, while the profile encloses a circular bubble of fluid in rigid body rotation at the ends of the gaps. Touching waves at the beginnings of the second gaps tend to the circular vortex wave on top of the limiting Crapper wave in the infinite vorticity limit, and the circular vortex wave on top of itself at the ends of the gaps. Touching waves for higher gaps accommodate more circular bubbles of fluid.

Equations

$${\mathbf{u}}_t+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\nabla P+(0,-g)\text{and}\nabla \cdot \mathbf{u}=0\text{in Ω(t),}$$

$$\omega =\nabla \times \mathbf{u}$$

$${\partial }_t+\mathbf{u}\cdot \nabla \text{ is tangent to ⋃tΓ(t)}\text{and}P=P_{atm}\text{ at Γ(t),}$$

$$\mathbf{u}(x,y,t)=(-\omega y,0)+\nabla \varphi (x,y,t)$$

$$\Delta \varphi =0\text{in Ω(t)}$$

$$\Omega (t)=\{(x,y)\in {\mathbb{R}}^2:y<0\},\mathbf{u}(x,y,t)=(-\omega y,0)\text{and}P(x,y,t)=P_{atm}-gy$$

$${\varphi }_t+\frac{1}{2}|\nabla \varphi {|}^2-\omega y{\varphi }_x+\omega \psi +gy=b(t)\text{at Γ(t)}$$

$$\varphi ,\psi \to 0\text{as y→−∞}\text{uniformly for x∈R}.$$

$${\varphi }_y=0\text{at y=−h}\text{for some h>0}$$

$$z(u+iv,t)=(x+iy)(u+iv,t)$$

$$z(u,t)=u+(\mathcal{H}+i)y(u,t)\text{for u∈R}.$$

$$\Gamma (t)=\{(u+\mathcal{H}y(u,t),y(u,t)):u\in \mathbb{R}\}.$$

$$\mathcal{H}{e}^{iku}=-i\text{sgn}(k)e^{iku}\text{for k∈Z}$$

$$W(u+iv,t)=(\varphi +i\psi )(z(u+iv,t),t)\text{for u+iv∈C−}.$$

$$W(u,t)=(1-i\mathcal{H})\varphi (u,t)\text{for u∈R}.$$

$$(1+\mathcal{H}{y}_u)y_t-y_u\mathcal{H}{y}_t-{\mathcal{H}}_d{\varphi }_u-\omega y{y}_u=0.$$

$$((1+\mathcal{H}{y}_u{)}^2+y_u^2)({\varphi }_t-\omega \mathcal{H}\varphi +gy-b(t))$$

$$-((1+\mathcal{H}{y}_u)\mathcal{H}{y}_t+y_u{y}_t){\varphi }_u+(y_u\mathcal{H}{y}_t-(1+\mathcal{H}{y}_u)y_t)\mathcal{H}{\varphi }_u$$

$$+\frac{1}{2}({\varphi }_u^2+(\mathcal{H}{\varphi }_u{)}^2)-\omega y((1+\mathcal{H}{y}_u){\varphi }_u-y_u\mathcal{H}{\varphi }_u)=0.$$

$$u\mapsto (u+\mathcal{H}y(u,t),y(u,t)),u\in \mathbb{R},\text{is injective}$$

$$((1+\mathcal{H}{y}_u{)}^2+y_u^2)(u,t)\ne 0\text{for any u∈R}.$$

$$⟨y(1+\mathcal{H}{y}_u)⟩=0$$

$$z_t=$$

$$W_t=$$

$$\alpha =\mathcal{P}\left(\frac{W_u-W_u^{\ast }+\frac{1}{2}\omega (z-z^{\ast })(z_u-z_u^{\ast })}{|z_u{|}^2}\right)\text{and}\beta =\mathcal{P}\left(\frac{|W_u{|}^2}{|z_u{|}^2}\right),$$

$$R=\frac{1}{z_u}\text{and}V=i\frac{W_u}{z_u}.$$

$$R_t$$

$$V_t$$

$${\varphi }^{\prime }=\mathcal{H}(\omega y{y}^{\prime }+c{y}^{\prime }).$$

$$(c^2+2b-2gy)((1+\mathcal{H}{y}^{\prime }{)}^2+(y^{\prime }{)}^2)-(c+\omega y(1+\mathcal{H}{y}^{\prime })-\omega \mathcal{H}(y{y}^{\prime }){)}^2=0.$$

$$\begin{gathered} (c^2+2b-2gy-{\omega }^2{y}^2)((1+\mathcal{H}{y}^{\prime }{)}^2+(y^{\prime }{)}^2) \\ -2\omega y(c-\omega \mathcal{H}(y{y}^{\prime }))(1+\mathcal{H}{y}^{\prime })=(c-\omega \mathcal{H}(y{y}^{\prime }){)}^2-{\omega }^2{y}^2(y^{\prime }{)}^2. \end{gathered}$$

$$(c^2+2b-2gy-{\omega }^2{y}^2)((1+\mathcal{H}{y}^{\prime }{)}^2+(y^{\prime }{)}^2)-2\omega y(c-\omega \mathcal{H}(y{y}^{\prime }))(1+\mathcal{H}{y}^{\prime }+i{y}^{\prime }).$$

$$(c^2+2b-2gy-{\omega }^2{y}^2)(1+\mathcal{H}{y}^{\prime }-i{y}^{\prime })-2\omega y(c-\omega \mathcal{H}(y{y}^{\prime })).$$

Full text

For strong positive vorticity, a ‘fold’ appears in the wave speed versus amplitude plane, and a ‘gap’ as the vorticity strength increases, bounded by two touching waves, whose profile contacts with itself at the trough line, enclosing a bubble of air. More folds and gaps follow for stronger vorticity.

Touching waves at the beginnings of the lowest gaps tend to the limiting Crapper wave as the vorticity strength increases indefinitely, while the profile encloses a circular bubble of fluid in rigid body rotation at the ends of the gaps. Touching waves at the beginnings of the second gaps tend to the circular vortex wave on top of the limiting Crapper wave in the infinite vorticity limit, and the circular vortex wave on top of itself at the ends of the gaps. Touching waves for higher gaps accommodate more circular bubbles of fluid.

1 Introduction

Stokes (1847, 1880) made formal but far-reaching considerations about periodic waves at the surface of an incompressible inviscid fluid in two dimensions, subject to the force of gravity, propagating a long distance at a practically constant velocity without change of form. For instance, he observed that the crests become sharper and the troughs flatter as the amplitude increases, and that the wave of greatest height, or the extreme wave, is distinguished by a $120^{\circ}$ peaking at the crest. It would be impossible to give here a complete account of Stokes waves. We merely pause to remark that in the irrotational flow setting, notable recent advances were based on a formulation of the problem as a nonlinear pseudodifferential equation, due to Babenko (1987) and others. For instance, Dyachenko et al. (2016); Lushnikov (2016); Lushnikov et al. (2017) revealed the complex singularities of almost extreme waves in meticulous detail.

The irrotational flow assumption is well justified in some situations. But rotational effects are significant in many others. Examples include wind driven waves, waves in a shear flow, and waves near a ship or pier. Constant vorticity is of particular interest. For instance, when waves are short compared with the characteristic lengthscale of vorticity, it is the value of the vorticity at the fluid surface which influences the waves (see Teles da Silva & Peregrine, 1988, for instance). Moreover, constant vorticity is representative of a wide range of physical scenarios. Examples include tidal currents, the alternating lateral movement of water in conjunction with the rise and fall of the tide, for which positive and negative constant vorticity are appropriate for the ebb and flood. By the way, here we distinguish negative vorticity for waves propagating downstream, and positive vorticity for upstream.

Dyachenko & Hur (2019) modified the Babenko equation and supplemented it with a scalar constraint, to permit constant vorticity and finite depth, and numerically solved by means of the Newton-GMRES method. Here we manipulate a change of variables to eliminate the Bernoulli constant from the modified Babenko equation and, hence, the scalar constraint. The associated linearized operator becomes self-adjoint, whereby efficiently handled by means of the Conjugate Gradient method. As a result, we achieve a high numerical precision for large amplitude waves for unprecedentedly large values of constant vorticity.

Simmen & Saffman (1985) discovered, and Dyachenko & Hur (2019) rediscovered, that for a large value of positive constant vorticity, a fold appears in the wave speed versus amplitude plane, and a gap of unphysical solutions as the value of the vorticity increases, bounded by two touching waves, whose profile contacts with itself at the trough line, enclosing a bubble of air. Here we find more folds and gaps as the value of the vorticity increases further. See Figure 5, for instance, for five folds and five gaps!

Moreover, Dyachenko & Hur (2019) discovered that touching waves at the beginnings of the lowest gaps tend to the limiting Crapper wave (see Crapper, 1957, for instance) as the value of positive constant vorticity increases indefinitely, revealing a striking and surprising link between rotational and capillary effects, while the profile encloses a circular bubble of fluid in rigid body rotation at the ends of the gaps. Here we find the ‘circular vortex wave’ on top of the limiting Crapper wave at the beginning of the second gap in the infinite vorticity limit, and the circular vortex wave on top of itself at the end of the gap. We find more circular bubbles of fluid at the boundaries of higher gaps in the infinite vorticity limit. See Figure 9, for instance, for a profile nearly enclosing five circular bubbles of fluid!

2 Formulation

We consider a two dimensional, infinitely deep and constant vorticity flow of an incompressible inviscid fluid, subject to the force of gravity, and the wave motion at the interface between the fluid and air. We assume for simplicity the unit density. Suppose that in Cartesian coordinates, the $x$ axis points in the direction of wave propagation, and the $y$ axis vertically upward. Suppose that the fluid occupies a region $\varOmega(t)$ in the $(x,y)$ plane at time $t$, bounded above by a free surface $\varGamma(t)$.

Let $\boldsymbol{u}=\boldsymbol{u}(x,y,t)$ denote the velocity of the fluid at the point $(x,y)$ and time $t$, and $P=P(x,y,t)$ the pressure. They satisfy the Euler equations for an incompressible fluid:

[TABLE]

where $g$ is the constant of gravitationals acceleration. Throughout, we assume that the vorticity

[TABLE]

is constant. We assume that there is no motion in the air and we neglect the effects of surface tension. The kinematic and dynamic conditions state that

[TABLE]

respectively, where $P_{atm}$ is the constant atmospheric pressure. We assume without loss of generality that $\varOmega$ and $\boldsymbol{u}$, $P$ are $2\pi$ periodic in the $x$ variable.

Suppose that

[TABLE]

for $(x,y)\in\varOmega(t)$, whence

[TABLE]

by the latter equation of (1). We merely pause to remark that for any $\omega\in\mathbb{R}$,

[TABLE]

solve (1) and (2) at all times. Here we assume that some external effects produce vorticity and we restrict the attention to waves propagating in a perturbed rotational flow. For nonconstant vorticity, (3) and (4) are inapplicable.

Substituting (3) into the former equation of (1) and using the latter equation of (2), we make an explicit calculation to arrive at

[TABLE]

for some function $b$, where $\psi$ is a harmonic conjugate of $\phi$. By the way, since $\phi$ and $\psi$ are determined up to the addition by arbitrary functions of $t$, we may take without loss of generality that $b(t)=0$ at all times.

The boundary condition in the infinite depth states that

[TABLE]

In the finite depth, we replace (6) by

[TABLE]

and proceed in like manner (see Dyachenko & Hur, 2019, for instance, for details). By the way, the effects of finite depth change the amplitude of a Stokes wave and others, but they are insignificant otherwise (see Dyachenko & Hur, 2019, for instance, for details).

2.1 Reformulation in conformal coordinates

We proceed as in Dyachenko & Hur (2019) and reformulate (1), (2) and (6) more conveniently in conformal coordinates. In the irrotational flow setting, we refer the interested reader to Dyachenko et al. (1996); Zakharov et al. (2002), among others. In what follows, we identify $\mathbb{R}^{2}$ with $\mathbb{C}$ whenever it is convenient to do so.

Suppose that

[TABLE]

conformally maps $\mathbb{C}^{-}:=\{u+iv\in\mathbb{C}:v<0\}$ of $2\pi$ period in the $u$ variable, to $\varOmega(t)$ of $2\pi$ period in $x$, at time $t$, such that $z(u+iv)\to u+iv$ as $v\to-\infty$. Suppose that (7) extends to map $\{u+i0:u\in\mathbb{R}\}$ to $\varGamma(t)$.

Since $z-(u+iv)$ is holomorphic in $\mathbb{C}^{-}$ and it is $2\pi$ periodic in the $u$ variable by hypothesis, a ‘Titchmarsh theorem’ (see Titchmarsh, 1986, for instance) implies that

[TABLE]

Therefore,

[TABLE]

Here and in the sequel, $\mathcal{H}$ denotes the periodic Hilbert transform, defined by

[TABLE]

and extended by linearity and continuity.

Moreover, let

[TABLE]

Since $W$ is holomorphic in $\mathbb{C}^{-}$ by (4) it is $2\pi$ periodic in the $u$ variable by hypothesis,

[TABLE]

Substituting (8) and (9) into the former equation of (2), we make an explicit calculation to arrive at

[TABLE]

They make what Dyachenko & Hur (2019), for instance, derived. Conversely, a solution of (10) gives rise to a solution of (1), (2) and (6), provided that

[TABLE]

See Dyachenko & Hur (2019), for instance, for details. By the way, an integration of (10a) over one period reveals that we may take without loss of generality that

[TABLE]

at all times. Here and in the sequel, $\langle f\rangle$ denotes the average of a $2\pi$ periodic function $f$ over one period.

Moreover, we solve (10) for $y_{t}$ and $\phi_{t}$ and, hence, $z_{t}$ and $W_{t}$ in the complex form. The result becomes, after an explicit calculation,

[TABLE]

where

[TABLE]

and $\mathcal{P}=\frac{1}{2}(1-i\mathcal{H})$. Here and in the sequel, the asterisk denotes complex conjugation. See Dosaev et al. (2017), for instance, for details. The imaginary part of the former equation of (11) and the real part of the latter make what Dyachenko & Hur (2019), for instance, derived. In the irrotational flow setting, they become what Dyachenko et al. (1996), for instance, derived.

Furthermore, let

[TABLE]

We solve (11) for $R$ and $V$, to arrive at

[TABLE]

See Dosaev et al. (2017), for instance, for details. In the irrotational flow setting, (12) becomes the so-called Dyachenko equation (see Zakharov et al., 2002, for instance), which is particularly convenient for numerical solution.

2.2 The Stokes wave problem

We turn the attention to traveling waves of (11), for which $z(u,t)-u$ and $W(u,t)$ are functions of $u-ct$ for some $c>0$, the speed of wave propagation, and $b(t)$ is constant. Under the assumption, we will go to a moving coordinate frame, changing $u-ct$ to $u$, whereby $t$ completely disappears. The former equation of (11) becomes

[TABLE]

Here and in the sequel, the prime denotes ordinary differentiation in the $u$ variable. The latter equation of (11), after an explicit calculation, simplifies to

[TABLE]

See Dyachenko & Hur (2019), for instance, for details.

Reformulation as the modified Babenko equation

We proceed as in Dyachenko & Hur (2019) and reformulate (13) more conveniently as an equation of Babenko kind. This makes use of that $\mathcal{H}+i$ determines the boundary value of a holomorphic and $2\pi$ periodic function in $\mathbb{C}^{-}$ up to the addition by a real constant.

We begin by rearranging (13) as

[TABLE]

Notice that $(c-\omega(\mathcal{H}+i)(yy^{\prime}))^{2}$ is the boundary value of a holomorphic and $2\pi$ periodic function in $\mathbb{C}^{-}$, whence so is

[TABLE]

Notice that $1/(1+\mathcal{H}y^{\prime}+iy^{\prime})$ is the boundary value of the holomorphic and $2\pi$ periodic function $1/z_{u}$ in $\mathbb{C}^{-}$, whence so is

[TABLE]

Therefore,

[TABLE]

up to the addition by a real constant. See Dyachenko & Hur (2019), for instance, for details. An integration over one period reveals that

[TABLE]

Indeed, $\langle\mathcal{H}f^{\prime}\rangle=0$ for any function $f$ and

[TABLE]

We emphasize that (14) is made up of polynomials of $y$, $y^{\prime}$ and their Hilbert transforms.

But for any $\omega,c$ and $b\in\mathbb{R}$, the left side of (14) maps a $2\pi$ periodic function to a $2\pi$ periodic function of mean zero, whereas the left side of (13) maps a $2\pi$ periodic function to a $2\pi$ periodic function, not necessarily of mean zero. In order to reconcile loss of information from (13) to (14), Dyachenko & Hur (2019) (see also Constantin et al., 2016) supplemented (14) with

[TABLE]

But the linearization of (14)-(15) with respect to $y$ and $b$ is not self-adjoint, because the linear part of (14) in $y$ depends on $b$, among others, whence the Conjugate Gradient or Conjugate Residual method is inapplicable. Dyachenko & Hur (2019) employed the Generalized Minimal Residual method, instead, and achieved success. But it becomes impractical for large amplitude waves for large values of positive constant vorticity.

Here we manipulate a change of variables

[TABLE]

for some $\eta_{0}\in\mathbb{R}$, so that $b=0$ but the linear and nonlinear parts of (14) remain invariant otherwise. This works, provided that

[TABLE]

Moreover, we perform a change of variables $y\mapsto y+\eta$ for some $\eta\in\mathbb{R}$, so that the constant part of (14) vanishes. The result becomes, by abuse of notation,

[TABLE]

An integration over one period reveals that

[TABLE]

Conversely, a solution of (14)-(15) and, hence, (16) gives rise to a periodic traveling wave of (1), (2) and (6), provided that

[TABLE]

See Dyachenko & Hur (2019), for instance, for details. By the way, $u+(\mathcal{H}+i)y(u)$ makes the boundary value of a conformal mapping (see (7)).

In what follows, we assume that $y$ is even (see Hur, 2007, for instance, for arbitrary vorticity).

To recapitulate, the Stokes wave problem in a constant vorticity flow is to seek a $2\pi$ periodic and even function $y$, satisfying (17), which solves (16) for some $\omega,c\in\mathbb{R}$. To compare, Dyachenko & Hur (2019) solved (14)-(15) for $y$ and $b$. But the associated linearized operator is not self-adjoint.

In the irrotational flow setting, (16) becomes the well-known Babenko equation (see Babenko, 1987, for instance). One may regard (16) as the ‘modified Babenko equation,’ permitting constant vorticity.

We emphasize that the linearization of (16):

[TABLE]

where $[f_{1},f_{2}]=f_{1}\mathcal{H}f_{2}-f_{2}\mathcal{H}f_{1}$, is self-adjoint, whereby we may employ the Conjugate Gradient method for efficient numerical solution. We achieve a high numerical precision for large amplitude waves for unprecedentedly large values of constant vorticity.

3 Numerical Method

3.1 Auxiliary conformal coordinates

Dyachenko et al. (2016) and others discussed that a Stokes wave defines the boundary value of a conformal mapping $z=z(w)$ for $w:=u+iv\in\mathbb{C}^{-}$ (see (7)), and an analytical continuation of $z$ to $\mathbb{C}$ has complex singularities in $[iv_{c},+i\infty)$, where $iv_{c}$ is the closest singularity of $z$ to $\mathbb{C}^{-}$; moreover,

[TABLE]

In the irrotational flow setting, $v_{c}\to 0$ as the profile tend to the wave of greatest height, or the extreme wave, presenting challenges in numerical solution. For instance, Dyachenko et al. (2016) took up to $2^{27}\approx 1.3\times 10^{8}$ uniform grid points in the $u$ variable for $v_{c}=5.93824419\dots\times 10^{-7}$. Lushnikov et al. (2017) devised an auxiliary conformal transform, involving trigonometric functions, to reduce the number of (nonuniform) grid points to $\sqrt{2^{27}}\approx 4.2\times 10^{4}$. See also Dyachenko & Hur (2019) in the constant vorticity setting.

But for a large value of positive constant vorticity, for instance, in Figure 1 are the graphs of $x=x(u)$ for two Stokes waves, whose profiles are in the last column of Figure 9; $(x(u),y(u))$, $u\in[-\pi,\pi]$, makes a conformal parametrization of a profile. The right panel suggests that $v_{c}$ is of the order of $10^{-24}$, whence the numerical method of Dyachenko et al. (2016); Lushnikov (2016); Lushnikov et al. (2017) would converge too slowly to be efficient.

Here we exploit conformal mappings to vertically slit regions in $\mathbb{C}$, involving Jacobi elliptic function, for further improvements. We refer the interested reader to Hale & Tee (2009), for instance, for details. For $L>0$, to be determined, let

[TABLE]

where $\operatorname{am}$ denotes the Jacobi amplitude and $K$ is the complete elliptic integral of the first kind with the elliptic modulus $\sqrt{1-L^{2}}$ (see Abramowitz & Stegun, 1964, for instance). We extend (19) to conformally map $\{q\in\mathbb{C}:-d<\text{Im}\,q<0\}$ of $2\pi$ period, to $\mathbb{C}^{-}$ of $2\pi$ period in the $w$ variable, where

[TABLE]

and $K^{\prime}$ is the complete elliptic integral of the first kind with the modulus $L$.

An explicit calculation reveals that (19) maps the Hilbert transform for the strip of depth $d$ in the $q$ variable to the Hilbert transform in the $w$ variable, so that (16) remains invariant in the $q$ variable. See also Silantyev (2019) in the irrotational flow setting.

Moreover, (19) maps uniform grid points in the $q$ variable over one period, to highly nonuniform in the $u$ variable. We wish to choose a suitable value of $L$, so that the inverse of (19) maps $iv_{c}$, the closest singularity to $\mathbb{C}^{-}$ in the $u$ variable, to $\text{Im}\,q=d$. Since

[TABLE]

it will enlarge considerably the region of analyticity in the $q$ variable when $v_{c}$ is small. For instance, if the closest singularity to $\mathbb{C}^{-}$ is of the order of $10^{-24}$ in the $u$ variable, the corresponding singularity is about $0.086$ in the $q$ variable, and in practice, it necessitates no more than $4096$ uniform grid points in the $q$ variable to resolve a numerical solution in quadruple precision.

{remark*}

Numerical evaluation of Jacobi elliptic functions in (19) is prone to catastrophic loss of precision in finite digit arithmetic, and it becomes more dramatic when $L$ decreases. Here we employ Landen transformations to evaluate the Jacobi amplitude, and the GNU MPFR library in arbitrary precision floating point arithmetic. We then truncate the result to quadruple precision arithmetic. We employ the arithmetic–geometric mean method to compute complete elliptic integrals. We refer the interested reader to Press et al. (1992); Bulirsch (1965), for instance, for details.

{remark*}

We use (20) to approximate the optimal value of $L$, but the exact value of $iv_{c}$, the closest singularity to $\mathbb{C}^{-}$ in the $u$ variable, is by and large unknown. When $L$ is smaller than optimal then the inverse of (19) maps $iv_{c}$ to $\pm a+id$ for some $a\in\mathbb{R}$, and $d$ depends on $L$ logarithmically. But if $L$ is larger than optimal then $iv_{c}$ is no longer mapped to $\text{Im}\,q=d$, and the Fourier spectrum in the $q$ variable widens rapidly.

3.2 The Newton-CG method

We write (16) in the operator form as

[TABLE]

where $y=y(u(q))$ and $q$ is in (19), and we solve it iteratively using the Newton method. Let

[TABLE]

where $y^{(0)}$ is an initial guess, to be supplied, and $\delta y^{(n)}$ solves

[TABLE]

$\delta\mathcal{G}(y^{(n)},c,\omega)$ is the linearization of $\mathcal{G}(y,c,\omega)$ with respect to $y$ and evaluated at $y=y^{(n)}$ (see (18)).

We take uniform grid points in the $q$ variable over one period. We approximate $y^{(n)}$ by a discrete Fourier transform and numerically evaluate using a fast Fourier transform. We numerically approximate polynomials of $y^{(n)}$, $(y^{(n)})^{\prime}$ and their Hilbert transforms likewise and, hence, $\mathcal{G}(y^{(n)},c,\omega)$ and $\delta\mathcal{G}(y^{(n)},c,\omega)$.

We numerically solve (22) using the Conjugate Gradient (CG) method (see Yang, 2009, 2010, for instance), in the Fourier space for the $q$ variable. Here we choose the linear part of $\mathcal{G}(y^{(n)},c,\omega)$ as a preconditioner. See Dyachenko et al. (2016), for instance, for details. We merely pause to remark that $\delta\mathcal{G}(y^{(n)},c,\omega)$ may not be positive definite, but the CG method converges nonetheless. To compare, the linearization of (14)-(15) is not self-adjoint, whence the CG method is inapplicable. Dyachenko & Hur (2019) employed the Generalized Minimal Residual Method (see Saad, 2003, for instance), instead, and achieved some success. The CG method is more powerful than the Generalized Minimal Residual method for a self-adjoint problem, and we achieve a high numerical precision for large amplitude waves for unprecedentedly large values of constant vorticity.

It turns out that we must supply an initial guess sufficiently close to a true solution of (21), so that the Newton-CG method will converge. We begin by taking $\omega=0$ and a numerical solution of small amplitude (see Dyachenko et al., 2016, for instance), and we continue the solution along in $\omega$, and $c$, if necessary, by taking the prior convergent solution as the initial guess and solving (21), until we arrive at a solution of small amplitude for the desired $\omega$. See Dyachenko & Hur (2019), for instance, for details. We then fix $\omega$ and continue the numerical solution along in $c$.

But the results in Section 4 (see also Dyachenko & Hur, 2019) reveal that for $\omega$ positive and large, a gap of ‘unphysical’ solutions appears in the wave speed versus amplitude plane; moreover, the gap becomes abnormally large in size as $\omega$ increases, so that it is impractical or even numerically prohibited to continue along in $c$ to trace out the gap. For instance, Figure 5 indicates that for $\omega\geqslant 2.4$, the wave speed reaches the order of hundreds in the gap. We then take a smaller $\omega$ and a physical numerical solution and we continue along in $\omega$, and $c$, if necessary, until we arrive at a physical solution for the desired $\omega$. See Dyachenko & Hur (2019), for instance, for details.

The results herein are carried out in quadruple precision floating point arithmetic. We say that the Newton method converges if the residuals $\lesssim N^{1/2}10^{-28}$, where $N$ is the number of uniform grid points in the $q$ variable. We take $N$ sufficiently large, so that the Fourier coefficient $\lesssim N^{-1/2}10^{-26}$ at the wave number $N$, to guarantee that the error in approximating a numerical solution by a discrete Fourier transform is inconsequential. To compare, Dyachenko & Hur (2019) required the residuals less than $10^{-13}$ and the Fourier coefficients less thans $10^{-12}$.

4 Results

Here we assume without loss of generality that $c$ is positive, and we permit $\omega$ positive or negative. Indeed, if $y$ and $c$ solve (16) for some $\omega$ then so do $y$ and $-c$ for $-\omega$. We take for simplicity that $g=1$.

In what follows, $s$ denotes the steepness, the crest-to-trough wave height divided by the period $=2\pi$, and $c$ measures the wave speed, for which $\langle y(1+\mathcal{H}y^{\prime})\rangle=0$. This is what Simmen & Saffman (1985); Teles da Silva & Peregrine (1988); Dyachenko & Hur (2019) chose, among others.

4.1 Folds and gaps

For a large value of positive constant vorticity, Simmen & Saffman (1985) discovered, and Dyachenko & Hur (2019) rediscovered, that the steepness increases, decreases and increases during the continuation of the numerical solution of the Stokes wave problem. Namely, a fold appears in the wave speed versus steepness plane; moreover, as the value of the vorticity increases, unphysical solutions make a gap in the wave speed versus steepness plane.

For instance, for $\omega=2.0$, the upper left panel of Figure 2 collects the wave speed versus steepness from the continuation of the numerical solution of (16). The right and lower left panels display, in the $(x,y)$ plane, the wave profiles of eight solutions, labelled by $A$ through $H$ along the $c=c(s)$ curve, in the order in which they appear in the numerical continuation. Here the troughs are at $y=0$. This is in excellent agreement with (Simmen & Saffman, 1985, Figure 10 and Figure 11), qualitatively and quantitatively. A caveat is that the vorticity of Simmen & Saffman (1985) differs in sign. We refer the interested reader to Simmen & Saffman (1985); Dyachenko & Hur (2019), among others, for $\omega$ negative, or positive but small.

The upper left panel reveals a fold in the $c=c(s)$ curve, along which $s$ increases and decreases. For $s$ small, for instance, for wave $A$, the profile is the graph of a single valued function. But we observe that the profile becomes more rounded as $s$ increases along the fold, so that overhanging waves appear, whose profile is no longer the graph of a single valued function; moreover, a touching wave follows, whose profile contacts with itself at the trough line, enclosing a bubble of air. For instance, wave $B$ is an almost touching wave. By the way, steep and rounded waves have relevance to backwash at the beach (see Teles da Silva & Peregrine, 1988, for instance).

We report that a numerical solution turns unphysical past the touching wave because (17a) no longer holds true. For instance, for wave $C$, the profile intersects itself and the fluid region overlaps itself. We observe that the wave profile becomes less rounded as $s$ decreases along the fold, so that another touching wave appears; moreover, a numerical solution turns physical past the touching wave. For instance, wave $D$ is unphysical, wave $E$ is an almost touching wave, and wave $F$ makes a physical solution past the touching wave. Together, a gap exists in the $c=c(s)$ curve, which consists of unphysical numerical solutions and which is bounded by two touching waves.

We merely pause to remark that the numerical method of Simmen & Saffman (1985) and others diverges in a gap. In stark contrast, the numerical method herein converges throughout. By the way, we discontinue a unphysical numerical solution when $c$ becomes abnormally large. Instead, we take a smaller $\omega$ and a physical numerical solution and we continue the solution along in $\omega$, and $c$, if necessary. See Dyachenko & Hur (2019), for instance, for details.

The insets of the upper left panel suggest that $s$ increases monotonically past the end of the fold, while $c$ experiences oscillations, like in the irrotational flow setting (see Dyachenko et al., 2016; Lushnikov et al., 2017, for instance). The lower left panel suggests that overhanging waves disappear as $s$ increases past the end of the fold; moreover, the inset reveals that the crests become sharper, like when $\omega=0$. We claim that an extreme wave appears ultimately, which is distinguished by peaking at the crest. One may not continue the numerical solution past the extreme wave because (17b) no longer holds true. For instance, wave $H$ is an almost extreme wave, and we calculate $s=1.72849105$. It is noticeably higher than the celebrated wave of greatest height or the extreme wave in the irrotational setting (see Dyachenko et al., 2016, for instance). But the extreme wave here seems no longer the wave of greatest height. For instance, for wave $E$, we calculate $s=1.78406376$. By the way, Simmen & Saffman (1985) estimated that the extreme wave steepness $\approx 1.727$. The method herein offers a significant improvement.

For a larger value of the positive constant vorticity, we take matters further and discover a second fold and a second gap in the wave speed versus steepness plane!

For instance, for $\omega=2.4$, in the upper left panel osf Figure 3 is the $c=c(s)$ curve from the continuation of the numerical solution of (16), and in the right and the lower left panels are the wave profiles of ten solutions at the indicated points along the $c=c(s)$ curve. The upper left panel reveals a fold, in the range $s=0$ to wave $D$, and a gap, bounded by waves $B$ and $C$, like Figure 2. The profile of wave $A$ is the graph of a single valued function. Wave $B$ is an almost touching wave near the beginning of the gap, and overhanging waves appear somewhere waves $A$ to $B$. Wave $C$ is an almost touching wave near the end of the gap. A numerical solution in the gap becomes unphysical because (17a) no longer holds true. Wave $D$ is near the end of the fold.

Moreover, the upper left panel reveals that $s$ increases from waves $D$ to $F$, decreases from waves $F$ to $H$, and increases past wave $H$, past the end of the fold. Together, another fold exists in the $c=c(s)$ curve.

The inset of the upper left panel suggests that $s$ increases monotonically past the end of the second fold, while $c$ experiences oscillations, like Figure 2 past the end of the fold. The lower left panel suggests that overhanging waves disappear as $s$ increases past the end of the second fold; moreover, the inset reveals that the crests become sharper, like when $\omega=2.0$. Therefore, we claim that an extreme wave appears ultimately as $s$ increases past the end of the second fold, and one may not continue the numerical solution past the extreme wave, because (17b) no longer holds true. For instance, wave $J$ is an almost extreme wave.

For $\omega=2.5$, in Figure 4 are the $c=c(s)$ curve and a selection of wave profiles at the indicated points along the $c=c(s)$ curve. The upper left panel reveals the lowest fold, in the range $s=0$ to wave $D$, and a gap, bounded by waves $B$ to $C$, and the second fold, from waves $D$ to $F$, like Figure 3.

We observe that the wave profile becomes more rounded as $s$ increases along the second fold, so that a touching wave appears, like along the lowest fold. For instance, wave $E$ is an almost touching wave. We observe that a numerical solution turns unphysical past the touching wave. For instance, for wave $F$, the profile intersects itself and the fluid region overlaps itself. We observe that the wave profile becomes less rounded as $s$ decreases along the second fold, so that another touching wave appears; moreover, a numerical solution turns physical past the touching wave, like along the lowest fold. For instance, wave $G$ is an almost touching wave and wave $H$ makes a physical solution past the touching wave. Together, a second gap exists in the second fold, consisting of unphysical numerical solutions and bounded by two touching waves.

We merely pause to remark that the numerical method of Simmen & Saffman (1985) and others diverges in a gap. In order to seek physical numerical solutions past the gap, Simmen & Saffman (1985) took $\omega=0$ and an almost extreme wave, manipulating that (17b) no longer holds true, and they continued the numerical solution along in $\omega$, and $c$, if necessary, until they arrived at an almost extreme wave for the desired $\omega$. They then continued the numerical solution along in $c$ until they arrived at an almost touching wave, past which their numerical method would diverge. This works, provided that there is one gap, but it is incapable of finding a second gap. To the best of the authors’ knowledge, a second gap reported herein is novel! By the way, the numerical method of Dyachenko & Hur (2019) converges in a gap, but it would take too much time to accurately resolve a numerical solution along a second fold.

The upper left panel suggests that $s$ increases monotonically past the second fold, while $c$ experiences oscillations, like Figure 3; moreover an extreme wave appears ultimately as $s$ increases past the second fold. For instance, wsave $J$ is an almost extreme wave.

To recapitulate, Figure 5 includes the wave speed versus steepness for several values of positive constant vorticity. For zero vorticity, Longuet-Higgins & Fox (1978) and others predicted that the wave speed experiences infinitely many oscillations while the steepness increases monotonically toward the wave of greatest height. Numerical computations (see Dyachenko et al., 2016; Lushnikov et al., 2017, for instance, and references therein) bear it out. The inset of the left panel reproduces it. For negative constant vorticity, the crests become sharper and lower (see Simmen & Saffman, 1985; Teles da Silva & Peregrine, 1988; Dyachenko & Hur, 2019, for instance). For a large value of positive constant vorticity, on the other hand, we find that the lowest oscillation of the wave speed deforms into a fold and, as the value of the vorticity increases, part of the fold transforms into a gap. For instance, for $\omega=1.74$, the left panel displays a fold and a gap. We refer the interested reader to Simmen & Saffman (1985); Dyachenko & Hur (2019), for instance, for details. We observe that the lowest fold and the lowest gap become larger in size as $\omega$ increases. For instance, for $\omega=2.0$, the fold and the gap are larger than $\omega=1.74$.

Moreover, for a larger value of the vorticity, for instance, for $\omega=2.4$, the figure reveals that the second oscillation of the wave speed deforms into another fold. As the value of the vorticity increases, for instance, for $\omega=2.5$, the right panel reveals that part of the second fold transforms into another gap. We observe that the second fold and the second gap become larger in size as $\omega$ increases. Furthermore, we find that higher folds and higher gaps develop in like manner. For instance, for $\omega=4.0$, the right panel reveals five folds and five gaps!

We observe that the steepness increases monotonically past all folds, although the wave speed experiences oscillations; moreover, overhanging waves disappear as $s$ increases past all folds, and the crests become sharper, like in the irrotational setting. Therefore, we claim that an extreme wave appears ultimately, which is distinguished by peaking at the crest. For instance, Figure 6 includes almost extreme waves for three values of positive constant vorticity, and bears it out. We calculate $\omega=2.0$, $s=1.72914844$; $\omega=2.4$, $s=2.54357018$, and $\omega=2.5$, $s=2.75939536$. Moreover, the figure suggests the extreme wave exhibits a $120^{\circ}$ peaking at the crest for any constant vorticity.

4.2 Touching waves in the infinite vorticity limit

Interestingly, Dyachenko & Hur (2019) discovered that for a sufficiently large value of positive constant vorticity, the touching waves at the boundaries of the lowest gap have capillary effects.

The left panel of Figure 7 displays the profiles of almost touching waves near the beginnings of the lowest gaps for four values of positive constant vorticity. We calculate $\omega=2.0$, $s=1.09717902$; $\omega=2.5$, $s=0.93652573$; $\omega=4.0$, $s=0.80492141$, and $\omega=14.0$, $s=0.73575899$. Indeed, numerical computations (see Dyachenko & Hur, 2019, for instance) suggest that $s$ decreases monotonically toward $\approx 0.73$ as $\omega\to\infty$.

In the irrotational setting, Crapper (1957) produced an exact solution of the capillary wave problem (in the absence of gravitational effects). He deduced that the crests become more rounded as the amplitude increases, opposite to gravity waves, and that the wave of greatest height is distinguished by the profile of a touching wave, for which $s\approx 0.730$. Moreover, the left panel reveals that for instance, for $\omega=14.0$, the profile of an almost touching wave matches remarkably well the limiting Crapper wave. Therefore, we claim that touching waves at the beginnings of the lowest gaps tend to the limiting Crapper wave as the value of positive constant vorticity increases unboundedly. This reveals a striking and surprising link between rotational and capillary effects.

In the right panel of Figure 7 are the profiles of almost touching waves near the ends of the lowest gaps for four values of positive constant vorticity. We calculate $\omega=2.0$, $s=1.78406173$; $\omega=2.5$, $s=1.70366636$; $\omega=4.0$, $s=1.41081893$, and $\omega=14.0$, $s=1.05512149$. Indeed, numerical computations (see Dyachenko & Hur, 2019, for instance) suggest that $s\to 1$ as $\omega\to\infty$. Teles da Silva & Peregrine (1988) numerically computed ‘pure rotational waves’ in a constant vorticity flow, in the absence of gravitational effects, and argued that a limiting wave is distinguished by a circular configuration of fluid in rigid body rotation. Moreover, the right panel reveals that for instance, for $\omega=14.0$, the profile of an almost touching wave is nearly circular. Therefore, we claim that touching waves at the ends of the lowest gaps tend to the ‘circular vortex wave’ as $\omega\to\infty$.

Moreover, in the left panel of Figure 8 are the profiles of almost touching waves near the beginnings of the second gaps, and the right panel near the ends of the gaps, for three values of positive constant vorticity. They lead to that touching waves at the beginnings of the second gaps tend to the circular vortex wave on top of the limiting Crapper wave as the value of positive constant vorticity increases unboundedly, while the circular vortex wave on top of itself at the end of the gaps in the infinite vorticity limit.

We claim that touching waves at the boundaries of higher gaps accommodate more circular bubbles of fluid in like manner! For instance, Figure 9 displays almost touching waves near the beginnings and ends from the lowest to fifth gaps for two large values of positive constant vorticity, and suggests that touching waves at the beginnings of the $n$-th gaps tend to $n-1$ circular vortex waves on top of the limiting Crapper wave, while $n$ circular vortex waves at the ends of the $n$-th gaps as $\omega\to\infty$.

Acknowledgements

VMH is supported by the National Science Foundation under the Faculty Early Career Development (CAREER) Award DMS-1352597. SD is supported by the National Science Foundation under DMS-1716822.