1. Introduction
Much of the theory of free-surface water waves assumes that the flow is irrotational, an assumption well justified in some applications. Also, in the absence of initial vorticity, boundaries or external forces, vorticity will remain zero for all times. Excellent surveys of travelling water waves can be found in Toland (Reference Toland1996), Strauss (Reference Strauss2010), Vanden-Broeck (Reference Vanden-Broeck2010) and Haziot et al. (Reference Haziot, Hur, Strauss, Toland, Wahlén, Walsh and Wheeler2022). However, rotational effects play a crucial role in numerous physical situations. For example, in any region where wind is blowing, there is a surface drift of the water and wave characteristics such as maximum wave height are sensitive to the velocity in the wind-drift boundary layer. See, for instance, Teles da Silva & Peregrine (Reference Teles da Silva and Peregrine1988) for further discussion.
A special case of free-surface waves with vorticity was found by Gerstner (Reference Gerstner1809), who derived an exact solution for periodic travelling waves in infinite depth for some non-zero vorticity. This was later rediscovered independently by Froude (Reference Froude1862), Rankine (Reference Rankine1863) and Reech (Reference Reech1869). Notably, no such closed-form solution exists for zero vorticity. Many authors have also considered the case of constant vorticity. Among the most striking in constant-vorticity flows are profiles with multi-valued height and even profiles which intersect themselves tangentially above the trough, enclosing an air bubble. By contrast, in an irrotational flow, the wave profile must be the graph of a single-valued function (see Hur & Wheeler (Reference Hur and Wheeler2022) and references therein).
Non-zero vorticity introduces substantial challenges in the mathematical treatment of travelling water waves. The stream function is no longer harmonic, whence complex analysis techniques are not directly applicable. Also, the fluid surface is a priori unknown, and like in all free-boundary problems, one must fix the free boundary. Dubreil-Jacotin (Reference Dubreil-Jacotin1934) proposed a semi-hodograph transformation, which interchanges the vertical coordinate (independent variable) with the stream function (dependent variable), successfully proving the existence of small-amplitude solutions. Constantin & Strauss (Reference Constantin and Strauss2004) took matters further and established the existence of large-amplitude solutions. The semi-hodograph transformation has the advantage that it accommodates arbitrary vorticity distributions. Its disadvantage is that it cannot permit internal stagnation points, critical layers or overhanging profiles. Nonetheless, it has been widely employed in analytical studies (see, for instance, Strauss (Reference Strauss2010) and references therein) and numerical investigations (Dalrymple Reference Dalrymple1977; Thomas Reference Thomas1990; Ko & Strauss Reference Ko and Strauss2008a ,Reference Ko and Strauss b ; Amann & Kalimeris Reference Amann and Kalimeris2018). Another approach is the so-called flattening transformation, for which the vertical coordinate is scaled to a fixed height. It can successfully describe interior stagnation points and critical layers, but it remains incapable of treating overhanging profiles (Wahlén Reference Wahlén2009).
A great deal of research on travelling water waves focuses on constant vorticity. This is because of its analytical tractability: the fluid domain can be conformally mapped from a fixed one, where the governing equations can be formulated in terms of quantities at the fluid surface, similarly to zero vorticity. Numerical studies (Simmen & Saffman Reference Simmen and Saffman1985; Teles da Silva & Pullin & Grimshaw Reference Pullin and Grimshaw1988; Peregrine Reference Teles da Silva and Peregrine1988; Vanden-Broeck Reference Vanden-Broeck1994, Reference Vanden-Broeck1996; Ribeiro Jr. et al. Reference Ribeiro, Milewski and Nachbin2017; Dyachenko & Hur Reference Dyachenko and Hur2019b ,Reference Dyachenko and Hur a ; Guan Reference Guan2020) and rigorous analysis (Constantin, Strauss & Vărvărucă Reference Constantin, Strauss and Vărvărucă2016; Kozlov, Kuznetsov & Lokharu Reference Kozlov, Kuznetsov and Lokharu2020; Hur & Wheeler Reference Hur and Wheeler2020, Reference Hur and Wheeler2022) (see also Haziot et al. (Reference Haziot, Hur, Strauss, Toland, Wahlén, Walsh and Wheeler2022) and references therein), employing a conformal transformation, have revealed remarkable phenomena, including interior stagnation points, critical layers and overhanging profiles (the latter of these were proven by Constantin, Strauss & Vărvărucă (Reference Constantin, Strauss and Vărvărucă2021) to occur only for waves with adverse currents, that is, currents that are sheared in the opposite direction to that of wave propagation). These features are expected to persist for more general vorticity distributions, motivating the development of numerical methods capable of capturing such behaviours.
Here we employ a conformal transformation introduced and analytically examined by Wahlén & Weber (Reference Wahlén and Weber2024, Reference Wahlén and Weber2023), which accommodates interior stagnation points, critical layers and overhanging profiles for arbitrary vorticity distributions. Unlike Wahlén & Weber (Reference Wahlén and Weber2024), who reduced the problem to a non-local equation, we numerically solve (local) partial differential equations in a fixed rectangular domain using a finite-difference scheme. The accuracy of our numerical method is validated through comparisons with known solutions for zero and constant vorticity, where strong agreement is observed. Additionally, novel solutions are presented for affine vorticity functions and a two-layer constant-vorticity scenario.
To the best of the authors’ knowledge, this is the first numerical study of Stokes waves in an arbitrary vorticity flow, capturing the effects of critical layers and overhanging profiles. The solutions obtained herein can provide a basis for stability analysis and other applications. While the present approach accounts only for gravitational effects, it can be extended to incorporate additional physical effects such as density stratification, surface tension and hydroelasticity.
2. Formulation
We consider periodic travelling waves at the free surface of an incompressible inviscid fluid in two dimensions under the influence of gravity in a rotational flow. Suppose that in Cartesian coordinates, the
$x$
axis points in the direction of wave propagation and the
$y$
axis opposite to gravity. In a frame of reference moving with the wave, the fluid flow appears steady and occupies a domain bounded above by a free boundary
$y=\eta (x)$
and below by a fixed boundary
$y=0$
. Let
\begin{align} \Omega _{(x,y)}=\{(x,y)\in \mathbb{R}^2: -L/2\leqslant x\leqslant L/2, 0\lt y\lt \eta (x)\}, \end{align}
where
$L\gt 0$
is the wavelength, and
\begin{align} H=\int _{-L/2}^{L/2} \eta (x)\,{\rm d}x \end{align}
the mean fluid depth. We introduce a stream function
$\psi (x,y)$
such that the velocity of the fluid is given by
$(\psi _y,-\psi _x)$
, where subscripts denote partial differentiation. The vorticity satisfies
\begin{align} \text{(vorticity)}=-\nabla _{(x,y)}^2\psi (x,y)=\gamma (\psi (x,y)), \end{align}
where
$\gamma$
is an arbitrary function of one variable. Throughout we assume it is single-valued. The relative mass flux is defined as
\begin{align} Q=\int ^{\eta (x)}_0 \psi _y(x,y)\,{\rm d}y, \end{align}
which is independent of
$x$
. The Stokes wave problem then takes the form
\begin{align} &\nabla _{(x,y)}^2\psi =-\gamma (\psi ) \rm \,\,\,\,\,\,\,\,\,\,\text{in $\Omega _{(x,y)}$}, \end{align}
\begin{align} &|\nabla _{(x,y)}\psi |^2+2g\eta =B \rm \,\,\,\,\,\,\,\,\,\,\text{at $y=\eta (x)$}, \end{align}
\begin{align} &\psi =Q \rm \,\,\,\,\,\,\,\,\,\,\text{at $y=\eta (x)$}, \end{align}
\begin{align} &\psi =0 \rm \,\,\,\,\,\,\,\,\,\,\text{at $y=0$}, \end{align}
where
$g$
is the constant of gravitational acceleration and
$B$
the Bernoulli constant. Throughout we assume the wave profile is symmetric about
$x=0$
(see figure 1).
While (2.5) accounts only for gravitational effects, it can be easily extended to include additional physical effects such as capillarity and hydroelasticity.
Schematic of the fluid domain in a reference frame moving at the wave speed, where the free surface is given by
$y=\eta (x)$
and the effects of non-trivial vorticity are represented by the background shear flow, illustrated with horizontal arrows.
The conformal parametrisation of the fluid domain in the rectangular region
$\Omega _{(\alpha ,i\beta )}$
. The formulation consists of two field equations and two sets of boundary conditions for the two unknowns
$\psi$
and
$Y$
.
In what follows, we identify
$\mathbb{R}^2$
with
$\mathbb{C}$
whenever convenient. Following Wahlén & Weber (Reference Wahlén and Weber2024, Reference Wahlén and Weber2023), we introduce
\begin{equation} x=X(\alpha ,\beta )\quad \text{and}\quad y=Y(\alpha ,\beta ), \end{equation}
which conformally maps
\begin{align} \Omega _{(\alpha ,\beta )}=\{\alpha +i\beta \in \mathbb{C}: -\pi \leqslant \alpha \leqslant \pi , -d\lt \beta \lt 0 \} \end{align}
for some
$d\gt 0$
, the conformal depth, to the fluid domain
$\Omega _{(x,y)}$
in the physical variables, satisfying
$\nabla _{(\alpha ,\beta )}^2 X, \nabla _{(\alpha ,\beta )}^2 Y=0$
. Additionally, they satisfy the Cauchy–Riemann equations,
$X_\alpha = Y_\beta$
and
$X_\beta =-Y_\alpha$
. A straightforward calculation reveals
\begin{align} \psi _x = \frac {1}{J}(\psi _\alpha Y_\beta - \psi _\beta Y_\alpha ) \quad \text{and}\quad \psi _y = \frac {1}{J}(\psi _\alpha Y_\alpha + \psi _\beta Y_\beta ), \end{align}
where
$J= |\nabla _{(\alpha ,\beta )} Y|^2$
, whence
$\nabla _{(\alpha ,\beta )}^2 \psi = J \nabla _{(x,y)}^2 \psi$
. Furthermore,
\begin{align} |\nabla _{(x,y)}\psi |^2 = \frac {1}{J} \psi _\beta ^2 \quad \text{at $\beta =0$}. \end{align}
Substituting these into (2.5), and imposing periodicity in the
$\alpha$
-variable, we arrive at
\begin{align} &\nabla ^2 \psi = -J \gamma (\psi ) \rm \,\,\,\,\,\,\,\,\,\,{in } \,\Omega _{\xi }, \end{align}
\begin{align} &\psi =Q \rm \,\,\,\,\,\,\,\,\,\,text{at } \beta =0, \end{align}
\begin{align} &\psi =0 \rm \,\,\,\,\,\,\,\,\,\,\text{at } \beta =-d, \end{align}
\begin{align} &\psi _\alpha = 0 \rm \,\,\,\,\,\,\,\,\,\,\text{at } \alpha =\pm \pi , \end{align}
and
\begin{align} &\nabla ^2 Y = 0 \rm \,\,\,\,\,\,\,\,\,\,\text{in } \Omega _{\xi }, \end{align}
\begin{align} &\frac {1}{J} \psi _\beta ^2 + 2 gY = B \rm \,\,\,\,\,\,\,\,\,\,\text{at } \beta =0, \end{align}
\begin{align} &Y=0 \rm \,\,\,\,\,\,\,\,\,\,\text{at } \beta =-d, \end{align}
\begin{align} &Y_\alpha = 0 \rm \,\,\,\,\,\,\,\,\,\,\text{at } \alpha =\pm \pi , \end{align}
where
$\nabla ^2=\nabla ^2_{(\alpha ,\beta )}$
(see figure 2).
If both the wavelength
$L$
and the mean fluid depth
$H$
are fixed, the relative mass flux
$Q$
and the conformal depth
$d$
must be allowed to vary. If only one of
$L$
or
$H$
is prescribed, on the other hand,
$d$
can remain fixed while
$Q$
varies. See § 4 for further discussion of bifurcation parameters.
Unlike the semi-hodograph transformation (Dubreil-Jacotin Reference Dubreil-Jacotin1934; Constantin & Strauss Reference Constantin and Strauss2004), the conformal mapping does not require streamlines to be graph-like in the
$x$
variable, thereby allowing for overhanging profiles. Also, internal stagnation points do not result in singularities in the mapping and therefore are admissible. However, care must be taken when stagnation points occur at the boundary of the fluid domain.
An important limiting scenario for travelling water waves involves the formation of a stagnation point at the fluid surface, characterised by a
$120^\circ$
angle at the crest. We highlight here the interesting result that for irrotational waves, this is dynamically only an apparent stagnation point of the flow. Despite the fact that the velocity vanishes there, no particle can rest there for a positive time, as shown for periodic waves by Constantin (Reference Constantin2006) and for solitary waves by Constantin & Escher (Reference Constantin and Escher2007). Since (2.6) is conformal, singularities cannot arise in the interior of the fluid domain, but they may still develop at the boundary. Specifically, mapping a
$120^\circ$
angle in
$z:=x+{\rm i}y$
to a straight line in
$\zeta :=\alpha +{\rm i}\beta$
would exhibit
$z\sim \zeta ^{3/2}$
to leading order in the vicinity of a stagnation point at the fluid surface and, hence, the second-order derivatives
$Y_{\alpha \alpha }$
and
$Y_{\beta \beta }$
become singular near the stagnation point. Consequently, our numerical method fails to converge for near-limiting and limiting solutions with an angle at the wave crest. See § 5 for further discussion of numerical limitations. An interesting direction for future research is to improve our numerical strategy to effectively resolve such singular behaviour.
It is interesting how our formulation simplifies under some assumptions. When the flow is irrotational, (2.10) is trivially satisfied by
$\psi =({Q}/{d})\beta$
because a unique conformal mapping from the fluid domain to a fixed rectangular region is a complex potential
$\phi +{\rm i}\psi$
, where
$\phi$
is a velocity potential such that the velocity
$=(\phi _x,\phi _y)$
. Therefore, solving (2.11) suffices to determine the solution in the absence of vorticity.
For a rigid wall, on the other hand, (2.11b
) is replaced by
$Y=H$
and (2.11) is trivially satisfied by
$Y=H+({H}/{d})\beta$
. Therefore, it remains to solve the vorticity equation. This corresponds to waves propagating in a shear flow in a uniform fluid channel. A promising extension of our numerical approach is to incorporate baroclinic vorticity arising from density stratification, which would be to include the free-surface effects in internal wave models (see e.g. Stastna (Reference Stastna2022) and references therein).
3. Linear theory
Before addressing nonlinear solutions, we examine the linear theory, which is relevant to small-amplitude perturbations of a flat surface. The dispersion relation for free-surface water waves in non-constant vorticity flows has been investigated in several studies, including Ehrnström et al. (2011), Karageorgis (Reference Karageorgis2012) and Kozlov, Kuznetsov & Lokharu (Reference Kozlov, Kuznetsov and Lokharu2014). Here we present a summary of the theory for completeness and outline a numerical solution method. It is important to note that the present analysis concerns steady waves and does not address stability, which would require working with the time-dependent problem where the vorticity function cannot be prescribed. The stability of Stokes waves in non-constant vorticity flows has been previous studied, for instance, in Hur & Lin (Reference Hur and Lin2008).
Let
\begin{align} \psi (x,y) = \psi _0(y) + \epsilon \psi _1(x,y)+\cdots \quad \text{and}\quad \eta (x) = \epsilon \eta _1(x)+\cdots \end{align}
for
$|\epsilon |\ll 1$
, where
$\psi _0(y)$
represents the background shear flow and is independent of
$x$
. Substituting these into (2.5) and collecting terms at leading order, we gather
\begin{equation} \left \{ \begin{aligned} &\frac {{\rm d}^2\psi _0}{{\rm d}y^2} = -\gamma (\psi _0) \quad \text{for}\quad 0\lt y\lt H \\ &\psi _0(H) = Q\quad \text{and}\quad \psi _0(0) = 0. \end{aligned}\right . \end{equation}
For some choices of
$\gamma (\psi )$
, an analytical solution is possible. In general, however, a numerical approach, such as using a shooting method, is required. At the order of
$\epsilon$
, we gather
\begin{align} \left \{\begin{aligned} &\nabla ^2 \psi _1= -\gamma '(\psi _0)\psi _1 && \text{for $0\lt y\lt H$} \\ &\psi _0^{\prime}(H) {\psi _1}_y + \psi _0^{\prime}(H) \psi _0^{\prime\prime}(H) \eta _1 + g \eta _1 = B_1 && \text{at $y=H$}\\ &\psi _1 + \psi _0^{\prime}(H) \eta _1= 0 && \text{at $y=H$}\\ &\psi _1=0 && \text{at $y=0$} \end{aligned}\right . \end{align}
for some constant
$B_1$
, where primes denote ordinary differentiation. Seeking solutions of the form
\begin{align} \psi _1(x,y) = f_1(y) \cos (kx) \quad \text{and}\quad \eta _1(x) = A_1\cos (kx) \end{align}
for some wavenumber
$k$
, we arrive at
\begin{align} &\,\,\, f_1^{\prime\prime} + (\gamma '(\psi _0)- k^2)f_1= 0 \quad \text{for $0\lt y\lt H$}, \end{align}
\begin{align} &\psi _0^{\prime}(H) f^{'}_{1}(H) + (\psi _0^{\prime}(H) \psi _0^{\prime\prime}(H) + g) A_1 = 0, \end{align}
\begin{align} &\qquad\qquad f_1(H) + \psi _0^{\prime}(H) A_1= 0, \end{align}
\begin{align} &\qquad\qquad\qquad\quad f_1(0) =0. \end{align}
To determine the dispersion relation, we proceed as follows. Given
$Q$
, we solve (3.2) for
$\psi _0(y)$
using a numerical shooting method. Substituting
$\psi _0(y)$
into (3.5), we solve (3.5a
) for
$f_1(y)$
subject to the boundary conditions (3.5c
) and (3.5d
) using a shooting method. It remains to determine the value of
$Q$
such that (3.5b
) is satisfied. This can be done graphically by plotting the residual of (3.5b
) as a function of
$Q$
or alternatively by employing an iterative solver such as the Newton–Raphson method.
In some special cases, the dispersion relation can be derived analytically. For instance, in an irrotational flow where
$\gamma (\psi )=0$
,
\begin{align} \psi _0(y) = \frac {Q}{H}y \quad \text{and} \quad f_1(y)= -A_1\frac {Q}{H}\frac {\sinh (ky)}{\sinh (kH)}, \end{align}
which leads to the well-known dispersion relation
$c^2 = {g}/{k}\tanh (kH)$
. For constant vorticity,
$\gamma (\psi )=\gamma _0$
, a straightforward calculation reveals
\begin{align} \psi _0(y) = cy-\frac {\gamma _0}{2}y^2 \quad \text{and} \quad f_1(y)= -A_1\left (\frac {c}{H}-\gamma _0\right )\frac {\sinh (ky)}{H\sinh (kH)}, \end{align}
where
$c={Q}/{H} + {\gamma _0}/{2} H$
. The dispersion relation then follows:
\begin{align} (c-\gamma _0 H )^2 + \frac {\gamma _0}{k}\tanh (kH)(c-\gamma _0 H ) - \frac {g}{k}\tanh (kH) = 0, \end{align}
which is consistent with the result, for instance, in Teles da Silva & Peregrine (Reference Teles da Silva and Peregrine1988).
Another instance where an analytical solution is possible is an affine vorticity function, that is,
\begin{align} \gamma (\psi ) = a\psi +b \quad \text{for some $a\gt 0$}\quad \text{for some $b\in \mathbb{R}$}. \end{align}
We follow the approach of Ehrnström et al. (2011) but with a different scaling. Solving (3.2), we obtain
\begin{align} \psi _0(y) = C\cos (\sqrt {a}y) + S\sin (\sqrt {a}y) - \frac {b}{a}, \end{align}
where
\begin{align} C&= \frac {b}{a} \quad \text{and}\quad S= \frac {Q+\frac {b}{a}(1-\cos (\sqrt {a}H))}{\sin {\sqrt {a}H}}, \end{align}
provided
$a\neq (n\pi /H)^2$
for
$n\in \mathbb{N}$
. If
$a=(n\pi /H)^2$
, then
$Q=0$
must hold. For notational convenience, we define
\begin{align} C_H= \cos (\sqrt {a}H) \quad \text{and} \quad S_H= \sin (\sqrt {a}H). \end{align}
Solving (3.5a ) and (3.5d ), we find
\begin{align} f_1(y) = \begin{cases} F \sin (\delta y), & a\gt k^2 \\ F y, & a=k^2 \\ F \sinh (\delta y), &a\lt k^2 \end{cases} \end{align}
for some constant
$F$
, where
$\delta ^2=a-k^2$
. To simplify the algebra, we focus on the case
$a\gt k^2$
, noting that similar expressions hold for
$a\leqslant k^2$
. It follows from (3.5c
) that
\begin{align} F = \frac {\sqrt {a}(CS_H-SC_H) }{\sin (\delta H)} A_1. \end{align}
Noting
$f^{'}_{1}(H) = \delta F \cos (\delta H)$
, (3.5b
) can be written as a quadratic equation for
$S$
:
\begin{equation} \begin{aligned}& (\delta C_H \cot (\delta H) + \sqrt {a} S_H ) C_H S^2 - \left(2\delta S_HC_H \cot (\delta H) + \sqrt {a} \left(S_H^2-C_H^2\right) \right)CS \\ &\quad + \delta \cot (\delta H) S_H^2C^2- \sqrt {a} C_HS_H C^2- \frac {g}{a} =0. \end{aligned} \end{equation}
We denote the two roots of this quadratic equation by
$S^\pm$
. Recalling the latter equation of (3.11), we obtain the dispersion relation
\begin{equation} Q^\pm = S^\pm S_H - C(1-C_H). \end{equation}
The algebra simplifies considerably when
$b=0$
, and
\begin{align} Q^2 = \frac {gS_H^2}{aC_H(\delta C_H \cot (\delta H) + \sqrt {a}S_H )}. \end{align}
In the long-wave limit as
$k\rightarrow 0$
,
\begin{equation} \sqrt {a} \cot (\sqrt {a}H) S^2 - \sqrt {a} CS - \frac {g}{a} =0 . \end{equation}
The dispersion relation forms asymptotes when the coefficient of
$S^2$
in (3.15) vanishes, that is, when
\begin{equation} \sqrt {a-k^2}\cot \big(\sqrt {a-k^2}\big) = -\sqrt {a} \tan (\sqrt {a}H). \end{equation}
The number of solutions for
$k$
to this transcendental equation is zero for sufficiently small
$a$
, but it increases as
$a$
increases. The existence of long waves depends on whether (3.18) has two real roots, which occurs if
\begin{align} b^2 + 4 \sqrt {a} g \cot \left(\sqrt {a}H\right)\gt 0. \end{align}
See § 5 for more details of the dispersion relation.
4. Numerical method
We describe how we numerically solve (2.10)–(2.11) using a finite-difference scheme. Our approach is similar to previous studies that reformulate the problem through the semi-hodograph transformation (Dalrymple Reference Dalrymple1977; Thomas Reference Thomas1990; Ko & Strauss Reference Ko and Strauss2008a ,Reference Ko and Strauss b ; Amann & Kalimeris Reference Amann and Kalimeris2018).
Since solutions are assumed to be symmetric about
$\alpha =0$
, we restrict computations to the half-domain
$\Omega _{(\alpha ,\beta )}^+$
where
$\alpha \geqslant 0$
. The boundary condition at
$\alpha =0$
is given by
$Y_\alpha =0$
, the same as the condition at
$\alpha =\pm \pi$
, so that our numerical method is applicable to symmetric and non-symmetric solutions, although for non-symmetric solutions computations span a full wavelength in
$\Omega _{(\alpha ,\beta )}^+$
rather than a half-wavelength.
While the method is formulated for periodic waves in finite depth, solitary waves can be approximated by taking the wavelength
$L$
sufficiently large such that further increases in
$L$
do not result in significant changes in the solution. Similarly, waves in infinite depth can be approximated by selecting a sufficiently large depth
$H$
.
We discretise
$\Omega _{(\alpha ,\beta )}^+$
into
$M$
equally spaced points in
$\alpha$
and
$N$
equally spaced points in
$\beta$
as
\begin{align} \alpha _m =\pi \frac {m-1}{M-1}\quad \text{and}\quad \beta _n =-d \frac {N-n}{N-1},\quad m=1,\ldots ,M, \quad n=1,\ldots ,N. \end{align}
The values of
$M$
and
$N$
are chosen depending on the computed solutions. For instance, for solitary waves approximated by long periodic waves, larger values of
$M$
are preferable to accurately resolve horizontal variations, while periodic waves in infinite depth necessitate larger values of
$N$
to adequately capture the flow motion at greater depths. Throughout § 5, the values of
$M$
and
$N$
used in computations are provided. We define
\begin{align} \Delta \alpha =\frac {\pi }{M-1} \quad \text{and} \quad \Delta \beta = \frac {d}{N-1}. \end{align}
Let
$Y_{m,n}$
and
$\psi _{m,n}$
denote the values of
$Y$
and
$\psi$
at the mesh point
$(\alpha _m, \beta _n)$
, respectively. Ghost points are introduced at the left and right boundaries, where
$\alpha _0=-\Delta \alpha$
and
$\alpha _{M+1}=\pi +\Delta \alpha$
. Applying the boundary conditions (2.10d
) and (2.11d
) at
$(\alpha _1,\beta _n)$
and
$(\alpha _M,\beta _n)$
, and employing second-order central differences, we obtain
$Y_{0,j}=Y_{2,j}$
and
$Y_{M+1,j}=Y_{M-1,j}$
, and similarly for
$\psi$
. The second-order central-difference approximations at the left and right boundaries lead to
\begin{equation} Y_{\alpha \alpha }(\alpha _1,\beta _n) \approx \frac {2(Y_{2,n}-Y_{1,n})}{\Delta \alpha ^2} \quad \text{and}\quad Y_{\alpha \alpha }(\alpha _M,\beta _n) \approx \frac {2(Y_{M-1,n}-Y_{M,n})}{\Delta \alpha ^2}. \end{equation}
The field equations are solved at internal mesh points
$(\alpha _m,\beta _n)$
, where
$m=1,\ldots ,M$
and
$n=2,\ldots ,N-1$
, using second-order central differences to approximate the derivatives. The boundary conditions (2.10c
) and (2.11c
) are enforced at
$(\alpha _m,\beta _1)$
for
$m=1,\ldots ,M$
. The boundary conditions (2.10b
) and (2.11b
) are imposed at
$(\alpha _m,\beta _N)$
for
$m=1,\ldots ,M$
, where
$\psi _\beta$
and
$Y_\beta$
are approximated using the second-order backward Euler equation, and
$Y_\alpha$
, for
$m=2,\ldots ,M-1$
, are approximated using second-order central differences, and set to zero at
$m=1$
and
$M$
(see figure 3).
Schematic of the discretisation of
$\Omega _{(\alpha ,\beta )}^+$
. The field equations and the boundary conditions are applied at mesh points corresponding to the coloured boxes. Crosses are mesh points, while circles show ghost points, which combined with (2.10d
) and (2.11d
) give second-order approximations (4.3).
This involves
$2MN+1$
unknowns – the values of
$\psi$
and
$Y$
at each mesh point along with the Bernoulli constant
$B$
– while providing only
$2MN$
equations. Additional constraints must be imposed to ensure a unique solution. This can be achieved by fixing bifurcation parameters. One approach is to fix the wavelength
\begin{align} L = 2\int _{0}^{\pi } Y_\beta (\alpha ,\beta )\,{\rm d}\alpha \Big|_{\beta ={\rm const.}}, \end{align}
which introduces an additional equation. Alternatively, one can fix the fluid depth
\begin{equation} H=\int _{0}^{\pi } Y(\alpha ,0)\,{\rm d}\alpha \quad \text{or} \quad H_0=Y(\pi ). \end{equation}
We introduce
$H_0$
because it is computationally less expensive to impose than the non-local equation for
$H$
. Furthermore, for solitary waves,
$H_0$
exactly corresponds to the fluid depth. Fixing the dimension of the mapping space, we can fix only one of
$L$
and
$H_0$
, and treat
$Q$
as an unknown. Additionally, we choose a bifurcation parameter, such as the wave amplitude
\begin{equation} A = Y_{1,N} - Y_{M,N}. \end{equation}
This results in
$2MN+2$
unknowns –
$\psi _{m,n}$
,
$Y_{m,n}$
,
$B$
,
$Q$
– which are matched by
$2MN+2$
equations –
$2MN$
from the field equations and boundary conditions, (4.6), and one of (4.4) or (4.5). If both
$L$
and
$H$
(or
$H_0$
) are prescribed, the dimension of the mapping space must be allowed to vary, that is,
$d$
is an unknown, leading to
$2MN+3$
unknowns –
$\psi _{m,n}$
,
$Y_{m,n}$
,
$B$
,
$Q$
,
$d$
. If neither
$H$
nor
$L$
is fixed, then both
$Q$
and
$d$
are prescribed, and there are
$2MN+1$
unknowns. Throughout § 5, we specify which physical parameters are fixed for each solution.
The Newton–Raphson method is employed for numerical solution. Convergence is achieved when the
$L^\infty$
-norm of the residuals falls below the tolerance
\begin{equation} {\rm TOL} = \frac {10^{-13}}{\max (\Delta \alpha ^2,\Delta \beta ^2)}, \end{equation}
which accounts for rounding errors introduced by double-precision floating-point arithmetic. Convergence is typically reached within a few iterations. The largest computational bottleneck is the construction of the Jacobian matrix at each iteration. While the linear part of the Jacobian remains fixed for a given
$d$
, the nonlinear terms must be updated. Fortunately, the sparsity of the Jacobian matrix allows for efficient inversion. Computational speed can be further improved by performing multiple iterations with the same Jacobian once the residuals are sufficiently small.
5. Results
Given the wide range of possible vorticity functions, a comprehensive treatment of all solutions is infeasible. Instead, we focus on a selected set of vorticity functions, chosen based on their physical relevance and connections to theoretical studies of rotational flows.
5.1. Zero vorticity
We begin by validating our numerical method against known results for solitary waves in an irrotational flow. To this end, we compare solutions from our finite-difference scheme with those obtained using a high-order series truncation method, for instance, in Doak & Vanden-Broeck (Reference Doak and Vanden-Broeck2017). Taking
$g=1$
and
$H_0=1$
, the far-field fluid velocity
$=\lim _{x\rightarrow \infty }|\psi _y(x,y)|$
, which, independent of
$y$
, is given by the non-dimensional Froude number, denoted
$F$
. We choose the conformal depth
$d$
sufficiently small such that further decreases in
$d$
do not alter the computed solutions significantly.
It is well established that small-amplitude solitary waves bifurcate from
$F=1$
and the amplitude increases with
$F$
until it reaches a limiting configuration which is distinguished by a stagnation point at the crest possessing a
$120^{\circ }$
angle. This limiting scenario was conjectured by Stokes (Reference Stokes1847, Reference Stokes1880) and later proved in Amick, Fraenkel & Toland (Reference Amick, Fraenkel and Toland1982) (see also Toland Reference Toland1978).
Solitary wave profiles for zero vorticity computed using the finite-difference scheme (red circles) and a high-order method (black curves). Solutions for (a)
$A=0.5$
and (b)
$A=0.8$
, showing excellent agreement between the two methods.
$(c)$
A wave of almost greatest height, where the finite-difference scheme fails to accurately capture the singularity behaviour with significant unphysical oscillations.
Figure 4 presents three wave profiles. Figures 4
$(a)$
and 4
$(b)$
show solutions with amplitudes
$A=0.5$
and
$A=0.8$
, respectively, demonstrating excellent agreement between the result from our finite-difference scheme (red circles) and that from the series truncation method (solid black curves). However, figure 4
$(c)$
shows the poor performance of the finite-difference scheme near the wave of greatest height, where
$A\approx 0.8332$
. It does not adequately resolve the singular behaviour near the crest, and spurious oscillations appear in the approximations of the derivatives.
Figure 5 presents a log–log plot of the error in the computed Froude number
$F$
versus the number of mesh points for
$A=0.5$
(solid black curve) and
$A=0.8$
(solid blue curve). The ‘true’ value of
$F$
is taken from the result using the series truncation method. Setting
$M=5N$
, we find that as
$M$
and, hence,
$N$
increase, the convergence of the numerical result is approximately quadratic in
$1/M$
. This is expected because of the use of second-order differences throughout the discretisation.
5.2. Constant vorticity
To further assess the accuracy and capability of our numerical scheme, we consider constant vorticity, for which (2.10) is no longer trivially satisfied.
Log–log plot of the error in the computed Froude number as a function of mesh sizing. For all solutions,
$N=M/5$
. The value
$F_{{true}}$
is obtained using a series truncation method. The black and blue curves correspond to the solutions shown in figures 4
$(a)$
and 4
$(b)$
. The dashed lines indicate quadratic convergence.
Solitary wave for constant vorticity
$\gamma =5$
for
$g=1$
,
$H_0=1$
and
$A=0.5$
.
$(a)$
Streamlines and (b) only the fluid surface. The black curves represent the result from our numerical method, while the red circles are computed in Guan (Reference Guan2020). The mesh sizes are
$M=2000$
and
$N=400$
. The bold black curves denote the fluid bed and a streamline that forms a critical layer.
We begin by examining a solitary wave for
$\gamma =5$
. We take
$g=1$
,
$H_0=1$
and
$A=0.5$
. Figure 6 presents the wave profile computed using our finite-difference scheme (black curves) alongside the result from a different numerical approach (red circles) solving for the conformal mapping of the irrotational correction to a linear shear current (e.g. Guan Reference Guan2020). Figure 6
$(b)$
displays the fluid surface, confirming strong agreement between the two methods. Figure 6
$(a)$
illustrates the streamlines throughout the fluid, revealing a critical layer with two stagnation points along the fluid bottom (bold black curves). Notably, such a solution would not be possible to compute using the semi-hodograph transformation.
Periodic travelling waves in infinite depth for
$\gamma =1$
and
$g=0$
. The black curves represent Crapper’s exact solutions, taken from Hur & Wheeler (Reference Hur and Wheeler2020). The red circles denote solutions computed using our numerical scheme, taking
$L=2\pi$
and
$d=7.5$
. Not all mesh points are displayed. (a)
$A=2$
and (b)
$A=4$
(see (5.1)).
$(c,d)$
The wave of maximum amplitude, with
$(d)$
providing a close-up of a touching region. The
$y$
axis has been shifted such that the interface is at
$y=0$
when
$x=\pm \pi$
. For all numerical solutions,
$M=200$
and
$N=3000$
.
Also, we consider zero gravity and infinite depth. Crapper (Reference Crapper1957) derived an exact solution for capillary waves in an irrotational flow of infinite depth, and a remarkable recent result established the equivalence between Crapper’s wave profiles and those in constant-vorticity flows in infinite depth in the absence of surface or body forces. This was hypothesised based on numerical and asymptotic evidence (Dyachenko & Hur Reference Dyachenko and Hur2019a
; Hur & Vanden-Broeck Reference Hur and Vanden-Broeck2020) and rigorously proved in Hur & Wheeler (Reference Hur and Wheeler2020). For
$\gamma =1$
and
$L=2\pi$
(
$d=\infty$
and
$g=0$
), the free surface can be parametrised as
$(X(\alpha ),Y(\alpha ))$
,
$-\pi \leqslant \alpha \leqslant \pi$
, where
\begin{align} X(\alpha ) + {\rm i} Y(\alpha ) &= \alpha - \frac {4{\rm i}C}{C+{\rm e}^{{\rm i}\alpha }}, \end{align}
where
$C\in [0,C_{max}]$
serves as an amplitude parameter. When
$C=0$
, the wave has zero amplitude (
$A=0$
). When
$C=C_{max}\approx 0.4547$
, on the other hand, the wave attains its maximum amplitude (
$A_{max}\approx 4.5851$
), for which the fluid surface intersects itself tangentially above the trough, enclosing an air bubble.
In figure 7, we compare solutions computed using our finite-difference scheme with the exact solutions in (5.1). While we cannot take true infinite depth, we approximate this by fixing the wavelength as
$2\pi$
and choosing a sufficiently large
$d$
such that further increases in
$d$
do not affect the numerical solution. For figure 7,
$d=7.5$
. We take
$M=200$
and
$N=3000$
. The solid black curves represent the exact solutions, while red circles denote our numerical result. Note that not all mesh points are displayed. Figures 7
$(a)$
and 7
$(b)$
depict solutions for amplitudes
$A=2$
and
$A=4$
, respectively. Figures 7
$(c)$
and 7
$(d)$
illustrate the wave of maximum amplitude, with the latter providing a close-up of a touching region. In all cases, the numerical and analytical solutions exhibit excellent agreement, further validating the accuracy of our numerical scheme.
5.3. Affine vorticity function
We turn attention to non-constant vorticity functions, with a particular emphasis on solutions that have thus far been numerically unattainable. The parameter space here is vast because the vorticity function can be arbitrary. Motivated by prior studies, we focus on two special cases: affine vorticity functions and two-layer constant vorticity.
The dispersion relation for
$\gamma (\psi )=a\psi + b$
for (a–d) different values of
$a$
and
$b$
.
Solutions for
$\gamma (\psi )=50\psi$
and
$k=1$
:
$Q=0.0413$
$(a)$
and
$Q=0.0499$
$(b)$
. All curves represent streamlines. The bold black lines denote the rigid wall and the fluid surface, while the bold blue lines indicate streamlines with saddle points. The solution in
$(a)$
is close to linear and exhibits two critical layers with saddle points at
$x=0$
. The solution in
$(b)$
approaches a limiting configuration with stagnation points at the fluid surface. There are multiple critical layers inside the flow, revealing an intricate flow structure caused by wave resonances. The solutions have
$M=500$
and
$N=200$
.
We begin with affine vorticity of the form
\begin{align} \gamma (\psi )=a\psi +b \quad \text{for some constants $a$ and $b$}, \end{align}
which can admit an arbitrary number of critical layers (Ehrnström et al. 2012). More specifically, examining the solution
$\psi _0(y)$
of (3.2), we observe that the velocity of the background shear flow oscillates in
$y$
with frequency
$\sqrt {a}$
. When exploring nonlinear solutions, it is helpful to first examine the dispersion relation in the linear theory. Figure 8 presents the dispersion relations for different values of
$a$
and
$b$
. As
$a$
increases, the number of asymptotes in
$Q(k)$
increases, their locations determined by (3.19). As seen in figure 8(a,c) (see also (3.16)),
$b=0$
results in the two solution branches with opposite values of
$Q$
. Indeed, when
$b=0$
, the equations are invariant under
$\psi \mapsto -\psi$
and
$Q\mapsto -Q$
, and the corresponding nonlinear solutions enjoy the same symmetry.
The asymptotes in the dispersion relation and, hence, multiple values of
$k$
for the same value of
$Q$
, suggest wave resonance. This occurs asymptotically when
$Q(k)=Q(nk)$
for a longer-wavelength mode with wavenumber
$k$
and for a shorter
$nk$
mode for some integer
$n\geqslant 2$
. Mathematical studies of resonance for capillary–gravity waves in an irrotational flow date back to Wilton (Reference Wilton1915).
Throughout this subsection, we set the length and time scales such that
$g=1$
and
$H=1$
, unless stated otherwise.
We begin by taking
$\gamma (\psi ) =50 \psi$
and compute the solution branch for
$k=1$
. Figure 9 presents the result. Figure 9
$(a)$
shows a small-amplitude solution, for which there are two critical layers inside the fluid. As the amplitude increases, wave resonances occur, resulting in additional local maxima and minima in the surface displacement. These resonances lead to intricate flow behaviours, including the formation of nested critical layers, as observed in a large-amplitude solution in figure 9
$(b)$
.
Our numerical result suggests that the limiting solution likely exhibits a stagnation point at the wave crest. However, as discussed earlier, our numerical method does not adequately resolve stagnation points at the fluid surface because the conformal mapping, more specifically its second-order derivatives, would become singular there. On the other hand, we observe that the fluid velocity near the crest decreases as one moves along the solution branch, supporting our conjecture that the solution branch would terminate at a solution with a stagnation point at the crest.
Solutions for
$\gamma (\psi )=50\psi +50$
and
$k=1$
:
$Q=-0.4859$
$(a)$
and
$Q=-0.5537$
$(b)$
. The solution in
$(a)$
resembles Wilton ripples. The
$k=1$
solution branch connects to a solution for
$k=6$
, shown in
$(b)$
. The solutions have
$M=500$
and
$N=200$
.
This is not the only limiting scenario by which solution branches terminate. For instance, taking
$a=50$
,
$b=50$
and
$k=1$
, we compute the solution branch bifurcating from
$Q=-0.2498$
. Figure 10 shows two representative solutions along the branch. In figure 10
$(a)$
, wave resonances occur, resulting in multiple local maxima and minima at the fluid surface. As one follows the solution branch, it reaches the solution shown in figure 10
$(b)$
, which lies on a solution branch with six periods in the computational domain. In other words, the
$k=1$
solution branch reaches a bifurcation point along the
$k=6$
solution branch.
Solution branch for
$\gamma (\psi )=50\psi$
and
$k=8$
.
$(a)$
The branch plotted in the
$(\min{q}, \max{\eta })$
plane. Bifurcating from an undisturbed interface for
$\max \eta =1$
, the amplitude increases monotonically as the crest speed decreases. (b,c) The solution marked by the cross in
$(a)$
, representing the farthest point along the branch where solutions could be obtained.
$(b)$
The streamlines in physical space;
$(c)$
the fluid velocity along the free surface as a function of
$x$
. The solutions were computed with
$M=N=500$
.
For a given
$a$
and
$b$
, wave resonances can be avoided by choosing
$k$
sufficiently large. For instance, figure 11 presents solutions with
$a=50$
,
$b=0$
and
$k=8$
. As shown in figure 8
$(c)$
, this choice of
$k$
ensures there are no resonant wavenumbers
$k\gt 8$
for the same value of
$Q$
. Under this condition, the wave amplitude increases monotonically along the branch, and the limiting solution exhibits a stagnation point at the crest. The critical layers here exhibit a far less complicated behaviour, forming Kelvin’s cat-eye streamline pattern, similar to that observed for constant vorticity (Ribeiro Jr. et al. Reference Ribeiro, Milewski and Nachbin2017).
Let
$q(x,\eta (x))$
denote the speed of the fluid at point
$(x,\eta (x))$
. We observe that for all solutions on the non-resonant branch,
$\min{q}$
occurs at
$x=0$
. Figure 11
$(a)$
depicts the solution branch in (
$\min q$
,
$\max{\eta }$
) space, while figure 11
$(b)$
presents the wave profile in physical space for the solution at the farthest point along the branch. In figure 11
$(c)$
, we plot
$q$
as a function of
$x$
for the solution in figure 11
$(b)$
. As we attempt to further increase the wave amplitude, the near-singular behaviour at the crest results in spurious oscillations, like those observed in an irrotational flow in figure 4
$(c)$
.
The amplitude of the near-limiting wave in figure 11 is much smaller than that for zero vorticity for the same parameter values for
$g$
,
$H$
and
$k$
. For constant vorticity, numerical studies (Dyachenko, Hur & Silantyev Reference Dyachenko, Hur and Silantyev2023) suggest that the limiting wave amplitude for negative vorticity is smaller than that for zero vorticity, whereas the amplitude for positive constant vorticity is larger. Furthermore, for positive constant vorticity, the wave amplitude can grow significantly larger than the irrotational counterpart. For instance, there are solution branches which approach a limiting configuration consisting of an infinite series of vertically stacked fluid bubbles in rigid-body rotation (Vanden-Broeck Reference Vanden-Broeck1996; Dyachenko & Hur Reference Dyachenko and Hur2019a
). This raises an interesting question as to what limiting amplitudes are achievable for affine vorticity functions.
Solutions for
$\gamma (\psi )=5\psi +5$
and
$k=1$
.
$(a)$
The wave profiles of near-limiting solutions along the branches bifurcating from
$Q=-2.55$
(red) and
$Q=-1.47$
(blue).
$(b)$
The background flow velocity. The red solution attains a near-limiting wave for a very small amplitude, where the background flow velocity is close to zero. The blue solution reaches a much larger amplitude, corresponding to the background flow near its maximum velocity at the fluid surface. The solutions were computed with
$M=1000$
and
$N=500$
.
We can make some predictions about the limiting amplitude given the background shear flow about which small-amplitude waves bifurcate. To illustrate this, we compute two branches of solutions for
$a=5$
,
$b=5$
and
$k=1$
, and present our findings in figure 12. Non-trivial solutions are found to bifurcate from
$Q=-1.47$
and
$Q=-2.55$
. In figure 12
$(b)$
, we show the background flow velocity
$\psi _0^{\prime}(y)$
for
$Q=-1.47$
(blue) and
$Q=-2.55$
(red). In figure 12
$(a)$
, we present the near-limiting solutions along these branches. We observe that the red solution reaches a limiting solution for a significantly small amplitude, where the background flow velocity is nearly zero at the fluid surface, whereas the limiting amplitude for the blue solution is much larger, corresponding to the background flow velocity near its maximum at the fluid surface. Given the oscillatory nature of the background shear flow, our results suggest that larger wave amplitudes are achieved when the maximum velocity is located close to the fluid surface.
Solutions for
$\gamma (\psi )=50\psi +50$
and
$k=6$
.
$(a)$
The solution branch in
$(Q, \max {\eta })$
space. Starting from a linear wave with
$Q=-0.385$
, the wave amplitude increases as
$Q$
decreases. Further along the branch, minimal changes in the streamlines occur as
$Q$
decreases further.
$(b)$
The solution marked by the black cross in
$(a)$
. Streamlines are depicted with the streamlines containing stagnation points. The red circles indicate the free surface of the solution marked by the red cross (
$Q=-20$
) in
$(a)$
. The solutions were computed with
$M=750$
and
$N=500$
.
Interestingly, this is not the only limiting scenario observed for non-resonant waves. For instance, we take
$a=50$
,
$b=50$
and
$k=6$
, and examine the branch of solutions bifurcating from
$Q=-0.385$
. One such solution along this branch, shown in figure 13
$(b)$
, features six wave periods. As the amplitude increases, the branch exhibits an increase in the magnitude of
$Q$
with almost no change in the overall flow structure. We plot this branch in
$(Q, \max{\eta })$
space in figure 13
$(a)$
. Computations were carried out up to
$Q=-65$
with
$M=750$
and
$N=500$
, beyond which further computations became impractical, requiring ever smaller steps in numerical continuation to achieve convergence. Additionally, for large values of
$Q$
, we had to relax the numerical tolerance by replacing
$10^{-13}$
by
$10^{-11}$
in (4.7). We attribute this to an increase in the order of the derivatives of
$\alpha$
and
$\beta$
, with
$\nabla ^2 \psi$
reaching values of the order of
$10^{4}$
. In this sense, the relative error remains comparable to that of small-amplitude solutions. In figure 13
$(b)$
, we show the largest-amplitude solution computed, with
$Q=-65$
. The red circles show the fluid surface and streamlines with saddle points for the solution with
$Q=-25$
, illustrating that there is little difference in the fluid surface and flow between these two solutions.
We note that all the limiting behaviours identified through our numerical computations are consistent with the theoretical predictions in Wahlén & Weber (Reference Wahlén and Weber2023).
5.4. Two-layer constant vorticity
Last but not least, we examine two layers of constant vorticity. Numerical computations for such a vorticity distribution can be found in Ko & Strauss (Reference Ko and Strauss2008a
,Reference Ko and Strauss
b
). Figure 12 in Ko & Strauss (Reference Ko and Strauss2008a
) concerns (after adjusting for differences of the choice of where
$\psi =0$
) solutions for
\begin{align} \gamma (\psi ) &= \begin{cases} 10, &-2\lt \psi \lt -1 \\ -10, & -1\lt \psi \lt 0. \end{cases} \end{align}
This creates two layers of constant-vorticity flows, with
$\gamma =10$
in the upper layer and
$\gamma =-10$
in the lower layer, introducing a discontinuity in the vorticity and, hence, in
$\nabla ^2 \psi$
. Since our computations involve second derivatives of
$\psi$
using finite differences, we instead seek a smooth vorticity function
\begin{align} \gamma (\psi ) = 10 \tanh \left (-40\left (\psi -\frac {Q}{2}\right )\right ). \end{align}
Solutions for (5.4) for
$g=9.8$
,
$H=0.6$
and
$k=1$
:
$Q=-2$
$(a)$
and
$Q=-3.25$
$(c)$
. Black lines indicate streamlines, and the colourbar represents vorticity. (b,d) The horizontal velocity at
$x=0$
for the solutions in (a,c). In
$(b)$
, the horizontal velocity is plotted as a function of
$\psi$
for comparison with Ko & Strauss (Reference Ko and Strauss2008a
), and in
$(d)$
as a function of
$y$
to highlight the broad region of near-stagnation. The solutions were computed with
$M=750$
and
$N=500$
.
Figure 14
$(a)$
compares these two vorticity functions for
$Q=-2$
, illustrating that (5.4) provides a smooth approximation of (5.3). To compare with figure 12 in Ko & Strauss (Reference Ko and Strauss2008a
), we fix
$g=9.8$
,
$H=0.6$
and
$k=1$
. The dispersion relation is shown in figure 14
$(b)$
. At
$k=1$
, we compute solutions bifurcating from
$Q=-1.573$
. As the amplitude increases, the value of
$Q$
decreases, and figure 15
$(a)$
shows the solution obtained at
$Q=-2$
, in excellent agreement with the result in Ko & Strauss (Reference Ko and Strauss2008a
).
Ko & Strauss (Reference Ko and Strauss2008a
) employ the semi-hodograph transformation, computing no further along the solution branch once the flow is at near stagnation inside the flow domain, which results in an almost singular change of variables. Figure 15
$(b)$
demonstrates this, where the negative velocity at
$x=0$
is plotted as a function of
$\psi$
, reproducing the result of Ko & Strauss (Reference Ko and Strauss2008a
). They conjecture that a stagnation point appears in the interior of the flow at the crest line as the amplitude is further increased.
Our solution method, which accommodates internal stagnation, allows us to compute solutions beyond interior stagnation. Figure 15
$(c)$
presents a large-amplitude solution, with the colourbar referring to the vorticity distribution. Notably, the region of zero vorticity, occupying only a small region at
$x=\pm \pi$
, expands to occupy a significant portion of the core of the wave. Figure 15
$(d)$
plots
$-u(x=0)$
as a function of
$y$
, revealing that the velocity nearly vanishes in this region. Continuing further along the branch, the amplitude continues to grow, and the near-stagnant bubble likewise grows in area.
This raises an intriguing question: what happens in the limiting solution as the vorticity function approaches (5.3)? Does this zero-vorticity region limit to a bubble inside a critical layer? A dedicated study of the two-layer vorticity configuration would be needed to explore this possibility. We note here that previous works have explored the dispersion relation for the case of free-surface waves with layers of constant (but different) vorticity in the fluid (Constantin Reference Constantin2012; Henry Reference Henry2013; Martin Reference Martin2014).
Solutions for (5.5) for
$g=9.8$
,
$H=0.6$
and
$k=1$
, where
$\gamma _0=5$
$(a)$
,
$10$
$(b)$
and
$15$
$(c)$
.
$(d)$
The magnitude of the horizontal velocity at the crest line, where the black, blue and red curves correspond to the solutions in (a–c). The solutions were computed with
$M=750$
and
$N=500$
.
The limiting solution of this branch remains unknown. As the amplitude increases further, our numerical method fails to successfully resolve the flow near the crest due to the stretching of mesh points by the conformal mapping. This highlights another limitation of conformal mappings: depending on the wave configuration and, hence, fluid domain, a conformal mapping from a rectangle with discretisation with equally spaced mesh points can lead to poor resolution in some regions.
We compared our result with figure 16 in Ko & Strauss (Reference Ko and Strauss2008b
) but found that even with a large number of mesh points,
$M=1300$
and
$N=500$
, the spacing near the crest became too spread out to maintain accuracy. For some solutions, solving in the
$(x,\psi )$
plane may provide better resolution. This issue could potentially be addressed with adaptive meshes or alternative coordinate mapping. We leave this for future work.
Another interesting observation is that the near-stagnant region appears only when there is a sign change in vorticity. For instance, figure 16 presents three solutions for the vorticity function
\begin{equation} \gamma (\psi ) = 10 \tanh \left (-40\left (\psi -\frac {Q}{2}\right )\right ) + \gamma _0. \end{equation}
Figure 16(a)–16(c) correspond to
$\gamma _0=5$
,
$10$
and
$15$
, respectively. The near-stagnant region is present in figure 16
$(a)$
but absent in figure 16(b,c).
Finally, we revisit the zero-gravity constant-vorticity solution computed in § 5.2. However, instead of constant vorticity, we consider a two-layer distribution of vorticity, connecting an irrotational bottom layer with a upper layer of constant vorticity of value unity. For this end, we choose
\begin{align} \gamma (\psi ) &= \frac {1}{2} { \left ( 1+\tanh { \left ( 40{ \left ( \psi -\frac {Q}{2} \right ) } \right ) } \right ) }. \end{align}
Fixing
$k=1$
and
$d=7.5$
, we show a large-amplitude solution in figure 17. The free-surface shape remains largely unchanged from the case of constant vorticity, where the free surface overhangs and nearly self-intersects. In this sense, the constant-vorticity upper layer seems largely unaffected by a change in the flow behaviour of the lower layer.
Solution for (5.6) with
$g=0$
,
$d=7.5$
,
$k=1$
and
$Q=8.551$
. The lower boundary is given by
$y=0$
. Black lines indicate streamlines, the blue streamline has a stagnation point that meets at a saddle and the colourbar represents vorticity. The profile remains largely unchanged from the case of constant vorticity throughout the whole fluid (shown with red circles).
6. Conclusion and future work
We have presented a numerical discretisation of the formulation introduced and analytically examined by Wahlén & Weber (Reference Wahlén and Weber2023) for periodic travelling water waves under gravity for an arbitrary distribution. Employing a conformal mapping rather than a semi-hodograph transformation, our approach allows for internal stagnation points, critical layers and overhanging profiles, broadening the range of computable solutions compared with previous studies (Dalrymple Reference Dalrymple1977; Thomas Reference Thomas1990; Ko & Strauss Reference Ko and Strauss2008b ,Reference Ko and Strauss a ; Amann & Kalimeris Reference Amann and Kalimeris2018). We have tested our numerical scheme against previously computed solutions for zero and constant vorticity, finding strong agreement.
We have uncovered novel solutions for affine vorticity functions, revealing intricate internal flow configurations, including multiple critical layers, theoretically predicted by Ehrnström et al. (2011, 2012). Furthermore, we have computed solutions where a smooth vorticity transition separates two layers of constant vorticity. When a sign change occurs in the vorticity, we observe that large-amplitude solutions develop near-stagnant regions in the transitional layer.
Some drawbacks of our numerical approach have been identified. While conformal mapping from a rectangular auxiliary domain allows for internal stagnation points, it does not permit singularities at the boundary. This presents challenges when a solution branch tends to stagnation at the fluid surface, because the local flow behaves like that inside a
$120^\circ$
corner, thereby the mapping is singular there. Additionally, since the method requires solving the problem throughout the fluid domain, rather than only at the fluid surface, as is for zero and constant vorticity, it can become computationally impractical to resolve some large-amplitude solutions.
Future work will involve a more comprehensive exploration of the solution space, focusing on both physically relevant and mathematically interesting vorticity functions. Additionally, modifications to the numerical method could allow for singular vorticity distributions in the bulk of the fluid and singularities at the boundaries, utilising a function-splitting method developed by Woods (Reference Woods1953). Finally, an intriguing potential extension is to consider vorticity induced by density variations, that is, baroclinic vorticity, where the density is given as a function of the stream function (see e.g. Long Reference Long1953). The methodology developed herein can be adapted to study free-surface solitary waves in a stratified fluid, replacing the rigid-wall assumption commonly used in previous studies (see Stastna (Reference Stastna2022) and references therein).
Acknowledgements
The authors thank the anonymous referees for their helpful suggestions.
Funding
A.D. would like to acknowledge funding from EPSRC NFFDy Fellowships (EPSRC grants EP/X028607/1). V.M.H. acknowledges funding from US NSF DMS-2407358.
Data availability
The codes used to produce the solutions seen in this paper (written in MATLAB) are available at https://github.com/alexdoak. Solutions can be provided upon request to the corresponding author.
Declaration of interests
The authors report no conflict of interest.
Open access
