1 Introduction
Nonlinear gravity–capillary waves are of interest due to their complexity and applications to wind–ocean coupling. Since the pioneering work of Longuet-Higgins (Reference Longuet-Higgins1989), considerable effort has been put into the numerical study of these waves of solitary type. It is well known that in deep water there are two fundamental branches of symmetric two-dimensional gravity–capillary solitary waves, including one family of waves with a positive free-surface elevation at their centre – denoted waves of elevation – and the other family of waves with a negative free-surface elevation at their centre – denoted waves of depression. Both branches bifurcate from infinitesimal periodic waves at the minimum of the phase speed where the group velocity is equal to the phase velocity (Vanden-Broeck & Dias Reference Vanden-Broeck and Dias1992). At small amplitude, these waves are in the form of wavepackets, and their envelopes can be approximated by the soliton solutions of the focussing cubic nonlinear Schördinger (NLS) equation (Dias & Iooss Reference Dias and Iooss1993).
Figure 1. (a) Bifurcation curve of the depression branch parameterised by the value at the middle point of waves, together with the typical wave
$1a$
(
$\unicode[STIX]{x1D702}(0)=-0.2433$
and
$c=1.4$
), the overhanging structure
$1b$
(
$\unicode[STIX]{x1D702}(0)=-1.48$
and
$c=0.9320$
) and the unphysical profile
$1c$
(
$\unicode[STIX]{x1D702}(0)=-1.50$
and
$c=0.9187$
). (b) Monotonic speed–energy relation for depression solitary waves.
Figure 2. (a) ‘Snake-like’ bifurcation of the elevation branch, which is parametrised by the value of the middle point. The curve starts from the bifurcation point
$(\sqrt{2},0)$
and follows the path
$2a\rightarrow 2b\rightarrow 2c\rightarrow 2d\rightarrow 2e\rightarrow 2f\rightarrow 2g\rightarrow 2h\rightarrow 2i\rightarrow 2j$
. The sharp nature of the second turning point (between
$2d$
and
$2e$
) and near overlapping points
$2b$
and
$2g$
are respectively shown in more details. The typical profiles for
$c=1.4$
(
$2d$
and
$2e$
) and
$c=1.36$
(
$2c$
and
$2j$
) are also shown. (b) Speed–energy bifurcation diagram of elevation solitary waves. Waves according to the segment between the turning point
$2B$
and the stationary point
$2C$
(solid line) are stable, while waves corresponding to the dotted line are unstable.
The speed–amplitude and speed–energy curves are both monotonic for the branch of depression (see Wang, Vanden-Broeck & Milewski (Reference Wang, Vanden-Broeck and Milewski2014) or figure 1 in this paper). As the translating speed
$c$
decreases, the depression solitary waves, which consist of a single trough at the centre, become steeper and gradually turn into an overhanging structure (i.e. multivalued wave profiles). The limiting configuration occurs when the two sides of the trough meet, enclosing a pendant-shaped bubble. The depression solitary waves were found to be stable by Calvo & Akylas (Reference Calvo and Akylas2002) using a numerical spectral analysis, however the result only applies to the single-valued wave profiles. On the other hand, Saffman (Reference Saffman1985) found a necessary condition for stability exchange for periodic gravity waves subject to superharmonic perturbations. His results show that the stability exchange can only occur at limiting points (including stationary points where the derivative is equal to zero and turning points where the derivative becomes infinite) of the total energy. As pointed out in his original paper, the argument can be easily generalised to include the surface tension term without essential modifications. For solitary waves, all the perturbations have the same horizontal scale as the fundamental waves, and hence are superharmonic. Therefore we can draw the conclusion that the depression solitary waves with an overhanging structure are also stable due to the monotonicity of the speed–energy curve (see figure 1
b). In this paper, the numerical evidence that the overhanging structure survives the collisions of solitary waves supports this conclusion.
Unlike the branch of depression solitary waves, stability exchanges may occur on the branch of elevation since its bifurcation diagrams feature different characteristics. Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014) traced the elevation branch and showed that there are many limiting points on the speed–energy curve where stability exchanges may occur according to Saffman’s theorem. It was found by Calvo & Akylas (Reference Calvo and Akylas2002) that the elevation branch is unstable between the bifurcation point and the first turning point, but gains spectral stability beyond this point, where the elevation waves resemble two depression waves placed side by side. These stability characteristics were confirmed in direct numerical simulations of the full equations by Milewski, Vanden-Broeck & Wang (Reference Milewski, Vanden-Broeck and Wang2010). These results raise the question as to whether stability exchanges can occur at other limiting points. More interestingly, after passing through the fourth turning point, the solitary-wave solutions, which consist of two troughs separated by two small dimples, appear (see
$2j$
in figure 2
a). Apart from the local shape between two basic troughs, the structure is similar to the known stable elevation waves (compare
$2j$
and
$2c$
in figure 2
a). One then can ask whether this structure is also stable. For the reasons given above, the stability of the elevation branch merits a revisit.
The existence of asymmetric solitary waves was originally predicted by Zufiria (Reference Zufiria1987) using the fifth-order Korteweg–de Vries (KdV) equation, a reduced model for small-amplitude gravity–capillary waves of finite depth when the Bond number
$B=\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C}gh^{2}$
(
$\unicode[STIX]{x1D70E}$
is the surface tension coefficient,
$\unicode[STIX]{x1D70C}$
is the fluid density,
$g$
is the gravitational acceleration and
$h$
is the depth of the fluid) is close to one-third. It is only fairly recently that the computations of asymmetric gravity–capillary solitary waves in the full Euler equations were carried out in deep water (Wang et al.
Reference Wang, Vanden-Broeck and Milewski2014). Apparently, there is the question of stability of these solutions. As emphasised by Yang & Akylas (Reference Yang and Akylas1997) and Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), these asymmetric waves are characterised by a multipacket structure, namely, they are formed by two or more adjacent wavepackets. Therefore we stand a chance of finding stable asymmetric waves if the constituent parts are all stable. However for a branch arising from an asymmetric solitary wave with the profile being composed of a depression wave and a stable elevation one, we give a negative answer in the present paper based on the linear eigenvalue analysis and direct numerical simulations.
On the experimental side, the depression gravity–capillary solitary waves can be excited by applying a localised pressure distribution moving with a constant subcritical speed; once the pressure is removed, a free depression solitary wave is shed and then quickly dissipates due to the viscous effect (see Longuet-Higgins & Zhang (Reference Longuet-Higgins and Zhang1997), Diorio et al. (Reference Diorio, Cho, Duncan and Akylas2009, Reference Diorio, Cho, Duncan and Akylas2011) for more details). Two-dimensional gravity–capillary travelling-wave solutions under a fully localised, constant-moving forcing were first computed by Vanden-Broeck & Dias (Reference Vanden-Broeck and Dias1992) in deep water. The linear stability of Vanden-Broeck & Dias’ solutions was later studied by Calvo & Akylas (Reference Calvo and Akylas2002), whereas in shallow water (Bond number bigger than one-third), the stability of the forced solutions was investigated by Grimshaw, Maleewong & Asavanant (Reference Grimshaw, Maleewong and Asavanant2009). The numerical studies of shedding of gravity–capillary solitary waves by forcing were carried out in the shallow water limit by Milewski & Vanden-Broeck (Reference Milewski and Vanden-Broeck1999) and Cho & Akylas (Reference Cho and Akylas2009) via the forced fifth-order KdV equation under the inviscid and viscous assumptions respectively, and in deep water by Akers (Reference Akers2012) via an inviscid cubic-truncation model. In this paper, the shedding of solitary waves in deep water is investigated using the direct time-dependent simulations of the fully nonlinear equations in the inviscid setting. Our results show that, besides the depression solitary waves, the stable elevation ones can also be excited by forcing, and the way to achieve this is to move two fully localised, well-separated pressures on the free surface with the speed slightly below the phase speed minimum and switch off the pressures simultaneously after a period of time.
2 Formulation
In this section, we recall the two-dimensional potential flow equations for a body of water of infinite depth. The equations can be reduced to a closed system of one dimension in space using the time-dependent conformal mapping technique, that maps the fluid domain to the lower half-plane where Laplace’s equation can be solved analytically. Following the procedure described by Milewski & Tabak (Reference Milewski and Tabak1999), the system is further recast to a single complex evolution equation and then numerically solved via the method of integrating factors for time integration. The free solitary-wave solutions are approximated by long periodic waves with flat tails, and computed using the pseudo-spectral method described in Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010), Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014).
2.1 Governing equations
Let
$x$
denote the horizontal coordinate,
$y$
the vertical direction, with
$y=0$
at the undisturbed interface between water and air, and
$t$
time. We consider a two-dimensional inviscid fluid of infinite depth flowing irrotationally, such that the velocity field can be expressed as the gradient of a potential function
$\unicode[STIX]{x1D719}$
. If the fluid is assumed to be incompressible, the equation that holds throughout the fluid is Laplace’s equation
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{xx}+\unicode[STIX]{x1D719}_{yy}=0,\quad \text{for }y<\unicode[STIX]{x1D702}(x,t), & & \displaystyle\end{eqnarray}$$
where the subscript denotes differentiation, and
$y=\unicode[STIX]{x1D702}(x,t)$
, which is also unknown, is designated as the free surface of the fluid. The kinematic and dynamic boundary conditions are to be satisfied on
$y=\unicode[STIX]{x1D702}(x,t)$
. The kinematic requirement that the free surface contains the same fluid particles is expressed by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{t}=\unicode[STIX]{x1D719}_{y}-\unicode[STIX]{x1D702}_{x}\unicode[STIX]{x1D719}_{x}. & & \displaystyle\end{eqnarray}$$
While the dynamic boundary condition reads
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{t}=-\frac{1}{2}\left(\unicode[STIX]{x1D719}_{x}^{2}+\unicode[STIX]{x1D719}_{y}^{2}\right)-g\unicode[STIX]{x1D702}+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x1D702}_{xx}}{(1+\unicode[STIX]{x1D702}_{x}^{2})^{3/2}}+\frac{P_{e}}{\unicode[STIX]{x1D70C}}, & & \displaystyle\end{eqnarray}$$
where
$g$
is gravitational acceleration,
$\unicode[STIX]{x1D70E}$
is surface tension,
$\unicode[STIX]{x1D70C}$
is fluid density and
$P_{e}$
represents the externally applied pressure at the free surface. Finally, the kinematic boundary condition at infinity is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{y}\rightarrow 0,\quad \text{as }y\rightarrow -\infty . & & \displaystyle\end{eqnarray}$$
The equations (2.1)–(2.4) form the complete system for the two-dimensional surface water-wave problem including both gravity and surface tension. Since gravity and surface tension are equally important in the present problem, one can non-dimensionalise the whole system by choosing
$$\begin{eqnarray}\displaystyle \left[\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}g}\right]^{1/2},\quad \left[\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}g^{3}}\right]^{1/4},\quad \left[\frac{\unicode[STIX]{x1D70E}^{3}}{\unicode[STIX]{x1D70C}^{3}g}\right]^{1/4}, & & \displaystyle\end{eqnarray}$$
as length, time and potential scales, respectively. Therefore the dynamic boundary condition (2.3) can be recast as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{t}=-\frac{1}{2}\left(\unicode[STIX]{x1D719}_{x}^{2}+\unicode[STIX]{x1D719}_{y}^{2}\right)-\unicode[STIX]{x1D702}+\frac{\unicode[STIX]{x1D702}_{xx}}{(1+\unicode[STIX]{x1D702}_{x}^{2})^{3/2}}+P_{e}. & & \displaystyle\end{eqnarray}$$
For a problem without external forcing, the total energy of the system is the sum of the kinetic and potential energies, which is given by
$$\begin{eqnarray}\displaystyle E=\frac{1}{2}\int _{-\infty }^{\infty }\,\text{d}x\int _{-\infty }^{\unicode[STIX]{x1D702}}\left(\unicode[STIX]{x1D719}_{x}^{2}+\unicode[STIX]{x1D719}_{y}^{2}\right)\,\text{d}y+\frac{1}{2}\int _{-\infty }^{\infty }\unicode[STIX]{x1D702}^{2}\,\text{d}x+\int _{-\infty }^{\infty }\left(\sqrt{1+\unicode[STIX]{x1D702}_{x}^{2}}-1\right)\text{d}x. & & \displaystyle\end{eqnarray}$$
Zakharov (Reference Zakharov1968) showed that the inviscid water-wave problem is a Hamiltonian system, and the correct Hamiltonian is simply the total energy given by (2.7) with surface variables
$\unicode[STIX]{x1D702}$
and
$\unicode[STIX]{x1D719}(x,\unicode[STIX]{x1D702},t)$
being canonical conjugates. In order to obtain the dispersion relation of gravity–capillary waves in deep water, we drop all the nonlinear terms in (2.2) and (2.6) and the external pressure in (2.6), and then solve the linearised system. It is easy to find the relation between the wave frequency
$\unicode[STIX]{x1D714}$
and the wavenumber
$k$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}^{2}=|k|\left(1+k^{2}\right),\quad c_{p}=\sqrt{\frac{1+k^{2}}{|k|}}, & & \displaystyle\end{eqnarray}$$
where
$c_{p}$
denotes the phase speed. It follows that the phase speed
$c_{p}$
attains its minimum
$\sqrt{2}$
at
$k=\pm 1$
. Vanden-Broeck & Dias (Reference Vanden-Broeck and Dias1992) showed that solitary waves bifurcate from the minimum of the phase speed and exist below this speed.
2.2 Numerical method
Based on the time-dependent conformal mapping technique (see Dyachenko, Lvov & Zakharov (Reference Dyachenko, Lvov and Zakharov1995), Dyachenko et al. (Reference Dyachenko, Kuznetsov, Spector and Sakharov1996), Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010) for details), the two-dimensional full Euler equations (2.1), (2.2), (2.4) and (2.6) can be reduced to a closed system of exact evolution equations. After mapping the physical fluid domain to the lower half-plane, the surface elevation
$\unicode[STIX]{x1D702}(\unicode[STIX]{x1D709},t)$
and the surface potential
$\unicode[STIX]{x1D719}(\unicode[STIX]{x1D709},t)$
can be parametrised by
$\unicode[STIX]{x1D709}$
, and the total energy in the transformed plane reads
$$\begin{eqnarray}\displaystyle E=\frac{1}{2}\int _{-\infty }^{\infty }\left[\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D702}^{2}+2\left(\sqrt{J}-x_{\unicode[STIX]{x1D709}}\right)\right]\text{d}\unicode[STIX]{x1D709}, & & \displaystyle\end{eqnarray}$$
where
$J\triangleq x_{\unicode[STIX]{x1D709}}^{2}+\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}^{2}$
is the Jacobian of the conformal map,
$x(\unicode[STIX]{x1D709},t)$
and
$\unicode[STIX]{x1D713}(\unicode[STIX]{x1D709},t)$
are associated with the harmonic conjugates of
$\unicode[STIX]{x1D702}(\unicode[STIX]{x1D709},t)$
and
$\unicode[STIX]{x1D719}(\unicode[STIX]{x1D709},t)$
respectively. More precisely,
$$\begin{eqnarray}\displaystyle x_{\unicode[STIX]{x1D709}}=1-{\mathcal{H}}\left[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}\right],\quad \unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}=-{\mathcal{H}}\left[\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}\right]. & & \displaystyle\end{eqnarray}$$
Here
${\mathcal{H}}$
represents the Hilbert transform defined by
$$\begin{eqnarray}\displaystyle {\mathcal{H}}\left[f\right](\unicode[STIX]{x1D709},t)=\frac{1}{\unicode[STIX]{x03C0}}\,\unicode[STIX]{x2A0D}_{-\infty }^{\infty }\,\frac{f(\unicode[STIX]{x1D709}^{\prime },t)}{\unicode[STIX]{x1D709}^{\prime }-\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}^{\prime }, & & \displaystyle\end{eqnarray}$$
where the integral is of the Cauchy principal form. The Fourier symbol of the Hilbert transform is
$\widehat{{\mathcal{H}}}=i\,\text{sgn}(k)$
where the hat denotes the Fourier transform and
$k$
is the wavenumber in Fourier space. Following Dyachenko et al. (Reference Dyachenko, Lvov and Zakharov1995), the Lagrangian function in the conformal plane takes the form
$$\begin{eqnarray}\displaystyle L & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D719}\left(\unicode[STIX]{x1D702}_{t}x_{\unicode[STIX]{x1D709}}-x_{t}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}\right)\,\text{d}\unicode[STIX]{x1D709}-E\nonumber\\ \displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\left[\unicode[STIX]{x1D702}_{t}\left(\unicode[STIX]{x1D719}x_{\unicode[STIX]{x1D709}}-{\mathcal{H}}\left[\unicode[STIX]{x1D719}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}\right]\right)-\frac{1}{2}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D713}-\frac{1}{2}\unicode[STIX]{x1D702}^{2}x_{\unicode[STIX]{x1D709}}-\left(\sqrt{J}-x_{\unicode[STIX]{x1D709}}\right)\right]\text{d}\unicode[STIX]{x1D709}.\end{eqnarray}$$
If we define
$\unicode[STIX]{x1D6E9}\triangleq \unicode[STIX]{x1D719}x_{\unicode[STIX]{x1D709}}-{\mathcal{H}}[\unicode[STIX]{x1D719}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}]$
,
$\unicode[STIX]{x1D6E9}$
and
$\unicode[STIX]{x1D702}$
are two canonical variables of the Hamiltonian system, and
$\unicode[STIX]{x1D719}$
can be recovered as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=\frac{\unicode[STIX]{x1D6E9}x_{\unicode[STIX]{x1D709}}+{\mathcal{H}}[\unicode[STIX]{x1D6E9}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}]}{J}. & & \displaystyle\end{eqnarray}$$
Then the evolution equations for
$\unicode[STIX]{x1D6E9}$
and
$\unicode[STIX]{x1D702}$
can be obtained by taking the variational derivatives of the total energy:
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x1D6FF}E}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6E9}},\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x1D6FF}E}{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D702}}. & & \displaystyle\end{eqnarray}$$
Based on the Hamiltonian formulation in the new coordinates, Saffman’s theorem of the stability exchange of water waves (Saffman Reference Saffman1985) can be extended to include multivalued solutions. However, the specific expressions of the Euler–Lagrangian equations (2.14) are rather cumbersome, therefore it makes sense to consider alternative variables to reduce the problem, even though these variables are not canonically conjugate. Following Dyachenko et al. (Reference Dyachenko, Kuznetsov, Spector and Sakharov1996) or Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010), we use
$\unicode[STIX]{x1D702}$
and
$\unicode[STIX]{x1D719}$
as dependent variables, then the dynamics of free surface in the transformed plane can be expressed as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{t}=\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}{\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]-x_{\unicode[STIX]{x1D709}}\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{t}=\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}^{2}-\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}^{2}}{2J}+\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}{\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]-\unicode[STIX]{x1D702}+\frac{x_{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}-\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}x_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}}{J^{3/2}}+P_{e}. & \displaystyle\end{eqnarray}$$
Following Milewski & Tabak (Reference Milewski and Tabak1999), we recast the system (2.15)–(2.16) as a single complex evolution equation. The numerical experiments are implemented on a periodic domain, where the Fourier transform is applicable in the spatial variable
$\unicode[STIX]{x1D709}$
, hence,
$$\begin{eqnarray}\displaystyle \left(\begin{array}{@{}c@{}}\widehat{\unicode[STIX]{x1D702}}\\ \widehat{\unicode[STIX]{x1D719}}\end{array}\right)_{t}+\left[\begin{array}{@{}cc@{}}0 & -|k|\\ 1+|k|^{2} & 0\end{array}\right]\left(\begin{array}{@{}c@{}}\widehat{\unicode[STIX]{x1D702}}\\ \widehat{\unicode[STIX]{x1D719}}\end{array}\right)=\left(\begin{array}{@{}c@{}}\widehat{{\mathcal{N}}}_{1}\\ \widehat{{\mathcal{N}}}_{2}\end{array}\right), & & \displaystyle\end{eqnarray}$$
with
${\mathcal{N}}_{1}$
and
${\mathcal{N}}_{2}$
being of the form
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{N}}_{1}=\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}{\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]+\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}\left(1-\frac{1}{J}\right)+\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}{\mathcal{H}}[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}]}{J}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{N}}_{2}=\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}^{2}-\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}^{2}}{2J}+\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}{\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]+\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}\left(\frac{1}{J^{3/2}}-1\right)+\frac{\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}{\mathcal{H}}[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}]-{\mathcal{H}}[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}]\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}}{J^{3/2}}+P_{e}. & \displaystyle\end{eqnarray}$$
Introducing the notations
$\widehat{p}=\widehat{\unicode[STIX]{x1D702}}+(i/c_{p})\widehat{\unicode[STIX]{x1D719}}$
and
$\widehat{q}=\widehat{\unicode[STIX]{x1D702}}-(i/c_{p})\widehat{\unicode[STIX]{x1D719}}$
, we can rewrite the system as
$$\begin{eqnarray}\displaystyle & \displaystyle \widehat{p}_{t}+i|\unicode[STIX]{x1D714}|\widehat{p}=\widehat{{\mathcal{N}}}_{1}+\frac{i}{c_{p}}\,\widehat{{\mathcal{N}}}_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \widehat{q}_{t}-i|\unicode[STIX]{x1D714}|\widehat{q}=\widehat{{\mathcal{N}}}_{1}-\frac{i}{c_{p}}\,\widehat{{\mathcal{N}}}_{2}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}$
and
$c_{p}$
are defined in (2.8). Using the fact that
$\unicode[STIX]{x1D702}$
and
$\unicode[STIX]{x1D719}$
are both real, these two equations are eventually equivalent, and
$\unicode[STIX]{x1D702}$
and
$\unicode[STIX]{x1D719}$
can be recovered from
$p$
alone with
$$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D702}}=\frac{1}{2}\;\left[\widehat{p}(k)+\widehat{p}(-k)^{\ast }\right],\quad \widehat{\unicode[STIX]{x1D719}}=\frac{c_{p}}{2i}\left[\widehat{p}(k)-\widehat{p}(-k)^{\ast }\right], & & \displaystyle\end{eqnarray}$$
where the asterisk indicates complex conjugation. Thus the problem is reduced to solving (2.20) which is a single complex evolution equation, and the method of integrating factors can be easily implemented for time integration (see Milewski & Tabak (Reference Milewski and Tabak1999) for more details).
In order to study the stability and dynamics of solitary waves, as a first step, we numerically seek fully localised travelling waves translating with speed
$c$
. We assume functions
$\unicode[STIX]{x1D702}$
and
$\unicode[STIX]{x1D719}$
to depend on
$\unicode[STIX]{x1D709}-ct$
, then it follows from (2.15) that
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\displaystyle -c\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}=\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}{\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]-x_{\unicode[STIX]{x1D709}}\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J},\\ \displaystyle -c{\mathcal{H}}\left[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}\right]=-x_{\unicode[STIX]{x1D709}}{\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]-\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
Solving the system (2.23) yields
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}=c\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}},\quad {\mathcal{H}}\left[\frac{\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D709}}}{J}\right]=c\left(\frac{x_{\unicode[STIX]{x1D709}}}{J}-1\right). & & \displaystyle\end{eqnarray}$$
Substituting (2.24) into the dynamic boundary condition (2.16) and noticing
$\unicode[STIX]{x1D719}_{t}=-c\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}$
, we obtain an integro-differential equation after some algebra
$$\begin{eqnarray}\displaystyle \frac{c^{2}}{2}\left(\frac{1}{J}-1\right)+\unicode[STIX]{x1D702}+\frac{\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}x_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}-x_{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}}{(J)^{3/2}}=P_{e}. & & \displaystyle\end{eqnarray}$$
It is noted that as
$P_{e}=0$
, the obtained solutions are free solitary waves. The nonlinear equation (2.25) was extensively used in computing gravity–capillary solitary waves (see, for example, Vanden-Broeck & Dias Reference Vanden-Broeck and Dias1992; Dias, Menasce & Vanden-Broeck Reference Dias, Menasce and Vanden-Broeck1996; Calvo & Akylas Reference Calvo and Akylas2002; Milewski et al.
Reference Milewski, Vanden-Broeck and Wang2010; Wang et al.
Reference Wang, Vanden-Broeck and Milewski2014). The total energy is a key parameter of free solitary waves since it is related to their stability characteristics according to Saffman’s theorem. In the transformed plane, the total energy takes the form
$$\begin{eqnarray}\displaystyle E=\frac{1}{2}\int _{-\infty }^{\infty }\,\left[c^{2}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}{\mathcal{H}}[\unicode[STIX]{x1D702}]+2\left(\sqrt{x_{\unicode[STIX]{x1D709}}^{2}+\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}^{2}}-x_{\unicode[STIX]{x1D709}}\right)+\unicode[STIX]{x1D702}^{2}x_{\unicode[STIX]{x1D709}}\right]\text{d}\unicode[STIX]{x1D709}. & & \displaystyle\end{eqnarray}$$
The solitary-wave solutions are approximated by long periodic waves with flat tails, therefore the Hilbert transform and derivatives can be efficiently computed in Fourier space using Fourier multipliers. We introduce a uniform mesh
$\unicode[STIX]{x1D709}_{i}=(i-1)(L/N)-(L/2)$
for
$i=1,2,\ldots ,N$
, and corresponding unknowns
$\unicode[STIX]{x1D702}(\unicode[STIX]{x1D709}_{i})$
, where
$L$
is the length of the domain and
$N$
is the number of grid points. The nonlinear integro-differential equation (2.25), which is discretised and evaluated in
$\unicode[STIX]{x1D709}-$
space, is solved via Newton’s method, and the solution is considered to have converged when the
$l^{\infty }$
-norm of the residual error is less than
$10^{-11}$
. Once one solitary wave is found, other solutions on the same branch can be computed via a numerical continuation method by choosing the wave speed or energy as a bifurcation parameter. This pseudo-spectral scheme was successfully implemented in Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010) for computing symmetric waves and in Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014) for finding asymmetric waves. Finally, the length of the computational domain
$L$
and the number of grid points
$N$
are both chosen sufficiently large so that, to the numerical accuracy we use, the solutions do not change when
$L$
and
$N$
are further increased.
3 Results
3.1 Solitary-wave dynamics
As shown in the aforementioned references, two basic families of symmetric gravity–capillary solitary waves, depression and elevation waves, are known to exist in deep water, and we briefly review some relevant results here. The bifurcation diagrams for depression and elevation branches together with typical wave profiles are presented in figures 1 and 2, from which we conclude that both branches bifurcate from infinitesimal periodic waves at the minimum of the phase speed
$c^{\ast }=\sqrt{2}$
, and exist at subcritical speeds. Figure 1 shows that both the amplitude (where ‘amplitude’ is defined as the centre elevation) and the energy (which is given by (2.26)) of depression waves are monotonic functions of the translating speed
$c$
. As
$c$
decreases the wave steepness increases, while the wave crest becomes more rounded and eventually the wave profile develops an overhanging structure at the trough (e.g. the profile
$1b$
). Further down the curve, the profile develops a point of contact with a ‘trapped bubble’ (which occurs at
$c\approx 0.9275$
) and even a self-intersecting structure (e.g. the profile
$1c$
). The solutions with self-intersections have no physical meaning but are mathematically admitted by the full Euler equations. In contrast to the branch of depression, the bifurcation of the elevation branch exhibits the characteristic ‘snake-like behaviour’: the speed–energy curve (figure 2
b) has stationary points where
$\text{d}E/\text{d}c=0$
(e.g.
$2C$
in figure 2
b), and turning points where
$\text{d}E/\text{d}c=\infty$
or equivalently
$\text{d}c/\text{d}E=0$
(e.g.
$2B$
in figure 2
b). Together with some wave profiles, the speed–amplitude curve of the elevation branch is shown in figure 2(a), which also features the snake-like behaviour. We stopped the computations after passing through the fourth turning point when two almost completely separated depression waves appear (see
$2j$
in figure 2
a). At this point, more accurate computations became prohibitive, presumably since the waves are far apart and their separation is sensitive to perturbations. We finally remark that there are infinite many branches of symmetric gravity–capillary solitary waves existing at finite amplitudes (the read is referred to the numerical computations by Dias et al. (Reference Dias, Menasce and Vanden-Broeck1996), Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), and the asymptotic analysis of Yang & Akylas (Reference Yang and Akylas1997)).
The stability characteristics of basic branches were examined by Calvo & Akylas (Reference Calvo and Akylas2002) via a linear spectral analysis and by Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2010) via direct numerical simulations. They found that all single-valued depression solitary waves are stable, while the elevation waves are unstable between the bifurcation point and the first turning point (
$2B$
in figure 2
b) but regain spectral stability beyond this point where the wave profiles feature two well-separated troughs (see
$2c$
in figure 2
a). The existence of stable solitary waves brings up the question of the excitation of these steady solutions. In laboratories, depression gravity–capillary solitary waves were generated under controlled conditions by Longuet-Higgins & Zhang (Reference Longuet-Higgins and Zhang1997), Diorio et al. (Reference Diorio, Cho, Duncan and Akylas2009, Reference Diorio, Cho, Duncan and Akylas2011), among others. The experimentalists excited the depression waves by blowing air towards the surface of the fluid and moving the air source with a speed close to the phase speed minimum. We present here an inviscid numerical analogue of these experiments. The equation (2.20) was simulated with a fully localised pressure distribution
$$\begin{eqnarray}\displaystyle P_{e}(x-ct)=-0.02\,\text{e}^{-(x-ct)^{2}},\quad \text{with }c=1.41. & & \displaystyle\end{eqnarray}$$
The time evolution of the transient response from a quiescent state, together with wave profiles, is illustrated in figure 3. Moving the frame of reference with the speed of the forcing, we observe a periodic shedding of depression solitary waves generated downstream of the pressure distribution. A similar phenomenon occurs for pure gravity waves in shallow water, where the KdV-type solitons can be periodically generated by a moving surface pressure, and the reader is referred to the comprehensive review paper by Keller (Reference Keller1985). The difference between these two cases is that the depression gravity–capillary solitary waves trail behind the pressure distribution, while the solitons for pure gravity waves advance upstream of the disturbance in processions. Despite being different in behaviour, two phenomena share an important core feature. In the transcritical range (namely, the speed slightly below the phase speed minimum for gravity–capillary waves in deep water, or the speed slightly above the phase speed maximum for pure gravity waves in shallow water) the dispersive effect is weak. When a fully localised pressure distribution moves with the speed in this range, there is neither steady forced solution nor linear mechanism to radiate energy away. Therefore the energy accumulates within the forced region, as a consequence, the amplitude of the local wave grows. When the amplitude reaches a certain level, the nonlinear effects start to become efficient. The subtle balance between weak dispersion and weak nonlinearity results in the formation of solitary waves moving at a different velocity from that of the forcing.
Figure 3. (a) Shedding of depression solitary waves resulting from moving a fully localised forcing with the constant speed
$c=1.41$
. (a) A forced solution in the
$\unicode[STIX]{x1D709}-t$
plane clearly indicates that solitary waves are generated periodically to the left of the forcing. (b) Snapshots in the physical domain at times
$t=500,1000,1500,2000$
from top to bottom respectively. The frame is moving with the forcing and only a portion of the computational domain is shown.
We now consider the generation of stable elevation solitary waves featuring two big troughs connected by a small dimple. Motivated by the structure of these waves, the pressure applied to the rest free surface consists of two identical, properly separated Gaussian functions
$$\begin{eqnarray}\displaystyle P_{e}(x-ct)=-0.03\,\left[\text{e}^{-(x-ct)^{2}}+\text{e}^{-(x-ct-10)^{2}}\right],\quad \text{with }c=1.405. & & \displaystyle\end{eqnarray}$$
Snapshots in figure 4(a) show that this special pressure impinging on the fluid surface can produce a disturbance in the form of an elevation wave. The pressure was then released after a period of time (
$t=80$
). The bottom snapshot in 4(a) and the top snapshot in 4(b) are shortly before (
$t=75$
) and after (
$t=100$
) the pressure was switched off respectively. The bottom snapshot in 4(b) indicates that the elevation solitary wave maintains its structure for long time in spite of the interaction with the background radiation field induced by the initial generation process.
Figure 4. Excitation of a stable elevation solitary wave due to two fully localised, well-separated pressures moving at
$c=1.405$
. The pressures are switched on at
$t=0$
and switched off at
$t=80$
. The solution (solid line) together with the pressures (dotted line) are shown in (a) at
$t=15,45,75$
from top to bottom respectively. After the pressures are turned off, the solution at times
$t=100,500,1000$
is shown in (b).
The numerical computations for generating figures 3 and 4 were performed on a periodic domain
$\unicode[STIX]{x1D709}\in [-250,250)$
, with the time step
$\unicode[STIX]{x0394}t=2\times 10^{-4}$
,
$N=20\,000$
Fourier modes, dealiasing with a buffer of
$N$
modes and no filtering. The character of the stable solution is not sensitive to the radiated waves wrapping around the domain.
The second dynamical computational experiment presented is the stability test of the depression solitary waves with an overhanging structure. As noted earlier, the numerical spectral analysis by Calvo & Akylas (Reference Calvo and Akylas2002) shows that depression waves are longitudinally stable. Unfortunately, their formulation is applicable only to the single-valued surface profiles, therefore the stability of the overhanging waves remains unknown. However this problem can be easily solved if one notices that the speed–energy curve is monotonic for the branch of depression and all perturbations are superharmonic in a sense for solitary waves. That’s because Saffman’s theorem indicates that for superharmonic perturbations stability exchanges can only occur at stationary points or turning points of the total energy. To provide more evidence for this assertion, we collided two depression waves: the larger slower wave was characterised by an overhanging structure, while the smaller faster wave served as a perturbation. Snapshots of the head-on collision in figure 5 show that both depression waves survive the collision, although they lose part of the energy, which eventually turns into a small wavepacket intercepted on the right-hand side of the larger solitary wave. The striking phenomenon is the persistence of the overhanging structure, although its shape slightly opens due to the energy loss during the collision process. In the case of overtaking collisions, the behaviours are different from those of head-on collisions, since the interacting time is much longer so that strong nonlinear effects are more likely to take place. An example is shown in figure 6, where only the larger solitary wave remains after the collision. As the smaller faster depression wave approaches the overhanging wave, it breaks up by transferring some energy to the overhanging structure and radiating the remaining part. We should emphasise that although the overhanging structure survives the overtaking collision, the gain of energy during the interaction results in a self-intersecting structure, therefore the solution satisfies the full Euler equations but has no physical significance. For the experiments of collisions, the equation (2.20) was integrated on a periodic domain
$\unicode[STIX]{x1D709}\in [-100,100)$
with the time step
$\unicode[STIX]{x0394}t=10^{-5}$
and
$N=4000$
Fourier modes.
Figure 5. Head-on collision between a solitary wave with an overhanging structure (
$1b$
in figure 1
a) and a small-amplitude depression wave (
$\unicode[STIX]{x1D702}(0)=-0.208$
and
$c=1.405$
): (a)
$t=0$
, (b)
$t=40$
, (c)
$t=50$
and (d)
$t=70$
. The collision is inelastic, but both depression solitary waves survive the collision.
Figure 6. Overtaking collision between an overhanging solitary wave (
$\unicode[STIX]{x1D702}(0)=-1.4676$
and
$c=0.94$
) and a small-amplitude depression wave (
$\unicode[STIX]{x1D702}(0)=-0.208$
and
$c=1.405$
) is shown in a frame of reference moving at the speed of the larger wave. From (a) to (d):
$t=0,200,250,300$
. Only the larger wave survives the collision. Note that an unphysical wave profile shows up after the collision.
3.2 Linear spectral analysis
In this section, using numerical spectral analysis and direct numerical simulations, we revisit the branch of symmetric elevation solitary waves and investigate the stability of a typical branch of asymmetric solitary waves. Following Tanaka (Reference Tanaka1985) or Calvo & Akylas (Reference Calvo and Akylas2002), we first outline the formulation of the linear eigenvalue problem for gravity–capillary waves. We introduce a frame of reference moving to the right with speed
$c$
, namely, we define
$\widehat{x}=x-ct$
and
$\widehat{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D719}-c\widehat{x}$
. Dropping hats, the dynamic boundary condition (2.6) can be rewritten as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{t}=-\frac{1}{2}\left(\unicode[STIX]{x1D719}_{x}^{2}+\unicode[STIX]{x1D719}_{y}^{2}\right)-\unicode[STIX]{x1D702}+\frac{\unicode[STIX]{x1D702}_{xx}}{(1+\unicode[STIX]{x1D702}_{x}^{2})^{3/2}}+\frac{c^{2}}{2}, & & \displaystyle\end{eqnarray}$$
while the kinematic boundary condition remains the same. As a consequence,
$$\begin{eqnarray}\displaystyle \left(\unicode[STIX]{x1D719}_{x},\unicode[STIX]{x1D719}_{y}\right)\rightarrow (-c,0),\quad \text{as }\sqrt{x^{2}+y^{2}}\rightarrow \infty . & & \displaystyle\end{eqnarray}$$
The steady solutions to the system (2.1), (2.2), (3.3) and (3.4) can be obtained by using hodograph transformation which exchanges the dependent variables
$\unicode[STIX]{x1D719}$
and
$\unicode[STIX]{x1D713}$
and independent variables
$x$
and
$y$
, and the Cauchy integral formula which solves Laplace’s equation in an explicit manner. Following Vanden-Broeck & Dias (Reference Vanden-Broeck and Dias1992), the horizontal and the vertical velocities can be expressed in terms of the complex potential
$f=\unicode[STIX]{x1D719}+i\unicode[STIX]{x1D713}$
and
$z=x+iy$
by the following formula
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{x}-i\unicode[STIX]{x1D719}_{y}=\left(\frac{\text{d}z}{\text{d}f}\right)^{-1}=\frac{1}{x_{\unicode[STIX]{x1D719}}+iy_{\unicode[STIX]{x1D719}}}, & & \displaystyle\end{eqnarray}$$
whereas
$x+iy$
is an analytic function of
$f$
in the upper half-plane (
$\unicode[STIX]{x1D713}>0$
). Applying the Cauchy integral formula to the function
$x_{\unicode[STIX]{x1D719}}+iy_{\unicode[STIX]{x1D719}}+(1/c)$
yields the relation between
$x_{\unicode[STIX]{x1D719}}$
and
$y_{\unicode[STIX]{x1D719}}$
at
$\unicode[STIX]{x1D713}=0$
$$\begin{eqnarray}\displaystyle x_{\unicode[STIX]{x1D719}}=-\frac{1}{c}+\frac{1}{\unicode[STIX]{x03C0}}\,\unicode[STIX]{x2A0D}_{-\infty }^{\infty }\,\frac{\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D719}^{\prime }}(\unicode[STIX]{x1D719}^{\prime })}{\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D719}}\,\text{d}\unicode[STIX]{x1D719}^{\prime }=-\frac{1}{c}+{\mathcal{H}}[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D719}}]. & & \displaystyle\end{eqnarray}$$
We rescale the Cauchy integral formula by defining
$\unicode[STIX]{x1D719}\triangleq c\,\unicode[STIX]{x1D709}$
, then the expression (3.6) becomes
$$\begin{eqnarray}\displaystyle x_{\unicode[STIX]{x1D709}}=-1+\frac{1}{\unicode[STIX]{x03C0}}\,\unicode[STIX]{x2A0D}_{-\infty }^{\infty }\,\frac{\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}^{\prime }}(\unicode[STIX]{x1D709}^{\prime })}{\unicode[STIX]{x1D709}^{\prime }-\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}^{\prime }=-1+{\mathcal{H}}[\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}]. & & \displaystyle\end{eqnarray}$$
It is not difficult to verify that the kinematic boundary condition is satisfied automatically by substituting (3.5) into (2.2), while the dynamic boundary condition (3.3) can be rewritten as
$$\begin{eqnarray}\displaystyle \frac{c^{2}}{2}\left(\frac{1}{x_{\unicode[STIX]{x1D709}}^{2}+\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}^{2}}-1\right)+\unicode[STIX]{x1D702}-\frac{\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}x_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}-x_{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D709}}}{(x_{\unicode[STIX]{x1D709}}^{2}+\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D709}}^{2})^{3/2}}=0. & & \displaystyle\end{eqnarray}$$
Evidently, (3.7) and (3.8) are coincident with (2.10) and (2.25), and this fact allows us to use the obtained travelling-wave solutions. When a fully localised steady solution is given, we consider small perturbations of this solution. The solitary waves are supposed to be single valued, as well as the perturbations. Therefore we can express the free-surface displacement
$\unicode[STIX]{x1D702}$
and the velocity potential
$\unicode[STIX]{x1D719}$
as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}(x,t)=Y(x)+\widetilde{\unicode[STIX]{x1D702}}(x)\text{e}^{\unicode[STIX]{x1D706}t},\quad \unicode[STIX]{x1D719}(x,y,t)=\unicode[STIX]{x1D6F7}(x,y)+\widetilde{\unicode[STIX]{x1D719}}(x,y)\text{e}^{\unicode[STIX]{x1D706}t}, & & \displaystyle\end{eqnarray}$$
where
$Y(x)$
and
$\unicode[STIX]{x1D6F7}(x,y)$
denote the profile of the free surface and the velocity potential of an exact solitary-wave solution, while
$\widetilde{\unicode[STIX]{x1D702}}\text{e}^{\unicode[STIX]{x1D706}t}$
and
$\widetilde{\unicode[STIX]{x1D719}}\text{e}^{\unicode[STIX]{x1D706}t}$
are small time-dependent perturbations superposed on these quantities respectively. Linearising the free boundary conditions (2.2) and (3.3), we obtain the eigenvalue problem for
$\unicode[STIX]{x1D706}$
,
$\widetilde{\unicode[STIX]{x1D702}}$
and
$\widetilde{\unicode[STIX]{x1D719}}$
as follows:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}\widetilde{\unicode[STIX]{x1D702}} & = & \displaystyle \left(\unicode[STIX]{x1D6F7}_{yy}-Y_{x}\unicode[STIX]{x1D6F7}_{xy}\right)\widetilde{\unicode[STIX]{x1D702}}-\unicode[STIX]{x1D6F7}_{x}\widetilde{\unicode[STIX]{x1D702}}_{x}-Y_{x}\widetilde{\unicode[STIX]{x1D719}}_{x}+\widetilde{\unicode[STIX]{x1D719}}_{y},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}\widetilde{\unicode[STIX]{x1D719}} & = & \displaystyle -\left(\unicode[STIX]{x1D6F7}_{x}\unicode[STIX]{x2202}_{x}+\unicode[STIX]{x1D6F7}_{y}\unicode[STIX]{x2202}_{y}\right)\unicode[STIX]{x1D6F7}_{y}\widetilde{\unicode[STIX]{x1D702}}-\left(\unicode[STIX]{x1D6F7}_{x}\unicode[STIX]{x2202}_{x}+\unicode[STIX]{x1D6F7}_{y}\unicode[STIX]{x2202}_{y}\right)\widetilde{\unicode[STIX]{x1D719}}-\widetilde{\unicode[STIX]{x1D702}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\widetilde{\unicode[STIX]{x1D702}}_{xx}}{(1+Y_{x}^{2})^{3/2}}-\frac{3Y_{x}Y_{xx}\widetilde{\unicode[STIX]{x1D702}}_{x}}{(1+Y_{x}^{2})^{5/2}}.\end{eqnarray}$$
Following Tanaka (Reference Tanaka1985) or Calvo & Akylas (Reference Calvo and Akylas2002), we introduce the magnitude of the velocity on the free surface,
$q=\sqrt{\unicode[STIX]{x1D6F7}_{x}^{2}+\unicode[STIX]{x1D6F7}_{y}^{2}}$
, and the angle between the velocity and the
$x$
-axis,
$\unicode[STIX]{x1D703}=\arctan \;(\text{d}Y/\text{d}x)$
. The free surface can be parameterised by the arc length
$s$
with
$s=0$
being consistent with
$\unicode[STIX]{x1D709}=0$
in (3.8). After some algebra, we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}\widetilde{\unicode[STIX]{x1D702}}=q\frac{\text{d}\widetilde{\unicode[STIX]{x1D702}}}{\text{d}s}-\frac{1}{\cos \unicode[STIX]{x1D703}}\frac{\text{d}\widetilde{\unicode[STIX]{x1D713}}}{\text{d}s}+\frac{1}{\cos \unicode[STIX]{x1D703}}\frac{\text{d}\left(q\cos \unicode[STIX]{x1D703}\right)}{\text{d}s}\widetilde{\unicode[STIX]{x1D702}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}\widetilde{\unicode[STIX]{x1D719}}=q\frac{\text{d}\widetilde{\unicode[STIX]{x1D719}}}{\text{d}s}-\left[q\,\frac{\text{d}(q\sin \unicode[STIX]{x1D703})}{\text{d}s}+1\right]\widetilde{\unicode[STIX]{x1D702}}+\frac{1}{\cos \unicode[STIX]{x1D703}}\frac{\text{d}}{\text{d}s}\left(\cos ^{2}\unicode[STIX]{x1D703}\frac{\text{d}\widetilde{\unicode[STIX]{x1D702}}}{\text{d}s}\right), & \displaystyle\end{eqnarray}$$
where
$\widetilde{\unicode[STIX]{x1D713}}$
is the harmonic conjugate of
$\widetilde{\unicode[STIX]{x1D719}}$
. Noticing
$\text{d}/\text{d}s=-q\text{d}/\text{d}\unicode[STIX]{x1D6F7}$
and
$\text{d}/\text{d}\unicode[STIX]{x1D709}=c\text{d}/\text{d}\unicode[STIX]{x1D6F7}$
, one finally rewrites the linear eigenvalue problem as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}\widetilde{\unicode[STIX]{x1D702}}=-\frac{c}{X_{\unicode[STIX]{x1D709}}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D709}}\left(\frac{X_{\unicode[STIX]{x1D709}}}{J}\right)\widetilde{\unicode[STIX]{x1D702}}-\frac{c}{J}\frac{\text{d}\widetilde{\unicode[STIX]{x1D702}}}{\text{d}\unicode[STIX]{x1D709}}+\frac{1}{X_{\unicode[STIX]{x1D709}}}\frac{\text{d}{\mathcal{H}}[\widetilde{\unicode[STIX]{x1D719}}]}{\text{d}\unicode[STIX]{x1D709}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}\widetilde{\unicode[STIX]{x1D719}}=-\left[\frac{c^{2}}{J}\frac{\text{d}(Y_{\unicode[STIX]{x1D709}}/J)}{\text{d}\unicode[STIX]{x1D709}}+1\right]\widetilde{\unicode[STIX]{x1D702}}-\frac{c}{J}\frac{\text{d}\widetilde{\unicode[STIX]{x1D719}}}{\text{d}\unicode[STIX]{x1D709}}-\frac{1}{X_{\unicode[STIX]{x1D709}}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D709}}\left(\frac{X_{\unicode[STIX]{x1D709}}^{2}}{J^{3/2}}\frac{\text{d}\widetilde{\unicode[STIX]{x1D702}}}{\text{d}\unicode[STIX]{x1D709}}\right), & \displaystyle\end{eqnarray}$$
where
$X_{\unicode[STIX]{x1D709}}=-1+{\mathcal{H}}[Y_{\unicode[STIX]{x1D709}}]$
and
$J=X_{\unicode[STIX]{x1D709}}^{2}+Y_{\unicode[STIX]{x1D709}}^{2}$
. We remark that under a different scaling, a similar formulation for left-going waves was previously used to study the stability of solitary waves for pure gravity waves in shallow water by Tanaka (Reference Tanaka1985) and for gravity–capillary waves in deep water by Calvo & Akylas (Reference Calvo and Akylas2002). In terms of the left- and right-going waves, the formulations given in Tanaka (Reference Tanaka1985) and presented here are essentially the same.
For consistency, the eigenvalue problem (3.14)–(3.15) is discretised at the uniform grid points of a periodic domain with the derivative and the Hilbert transform being approximated by the spectral matrices. Specifically, the square matrix of the discrete Hilbert transform, denoted by
$\mathbb{H}$
, is given by
$$\begin{eqnarray}\displaystyle \mathbb{H}(i,j)=\left\{\begin{array}{@{}ll@{}}\displaystyle -\frac{2}{N}\,\cot \left(\frac{\unicode[STIX]{x03C0}(i-j)}{N}\right),\quad & |i-j|\not \equiv 0\;(\text{mod}\;2),\\ 0,\quad & |i-j|\equiv 0\;(\text{mod}\;2),\end{array}\right. & & \displaystyle\end{eqnarray}$$
and the first-derivative matrix
$\mathbb{D}$
is given by
$$\begin{eqnarray}\displaystyle & \mathbb{D}(i,j)=\left\{\begin{array}{@{}ll@{}}\displaystyle (-1)^{i-j}\,\frac{\unicode[STIX]{x03C0}}{L}\,\cot \left(\frac{\unicode[STIX]{x03C0}(i-j)}{N}\right), & |i-j|\not \equiv 0\;(\text{mod}\;N),\\ \displaystyle 0, & \displaystyle |i-j|\equiv 0\;(\text{mod}\;N),\end{array}\right. & \displaystyle\end{eqnarray}$$
where
$L$
is the length of the computational domain and
$N$
is the number of grid points (see Trefethen (Reference Trefethen2000) for detailed derivations). Therefore, the original eigenvalue problem is approximated by a matrix eigenvalue problem, which can be solved using the standard QR algorithm in MATLAB.
To validate our numerical code, we checked the results produced by our program against the results in Calvo & Akylas (Reference Calvo and Akylas2002). For the elevation gravity–capillary solitary waves, Calvo & Akylas (Reference Calvo and Akylas2002) considered the stability exchange near the first turning point, and found that the waves are unstable before that point. For
$\unicode[STIX]{x1D6FC}$
(which is defined as
$1/c^{4}$
) equals
$0.26$
,
$0.3$
and
$0.4$
, the growth rate of the most unstable mode is respectively
$0.009$
,
$0.0153$
and
$0.0077$
. While our numerical results read
$0.009$
,
$0.0154$
and
$0.01096$
after being converted to their scaling. The small discrepancies are most likely due to the different way to approximate infinitely long solitary waves, and the different matrices of the Hilbert transform and derivatives arising therefrom. The profile of the most unstable mode for
$\unicode[STIX]{x1D6FC}=0.3$
was also checked (see figure 7), which is in excellent agreement with Calvo & Akylas’ result (figure 10 in their paper), confirming the accuracy of our numerical method. Throughout this work, the linear eigenvalue problem (3.14)–(3.15) was solved with
$L=200$
and
$N=8000$
for all waves we tested.
Figure 7. The most unstable normal mode for the elevation wave with
$\unicode[STIX]{x1D6FC}=0.3$
(which corresponds to
$c=1.3512$
). The profiles are modified in scale in order to compare with figure 10 of Calvo & Akylas (Reference Calvo and Akylas2002).
Table 1. Instability growth rate
$\unicode[STIX]{x1D706}$
of elevation solitary waves for
$c=1.36$
and
$c=1.4$
. The value of the middle point and the energy of the wave are listed for reference.
We start to explore the stability characteristics of the elevation branch near other limiting points. In figure 2(b), away from the turning point
$2B$
no instability can be detected until the stationary point
$2C$
, beyond which the elevation waves become unstable again. Table 1 shows the instability growth rate of unstable elevation waves at
$c=1.36$
and
$c=1.4$
. The most striking result is the instability of the waves after the fourth turning point, whose profiles feature two big troughs connected by two small dimples (see
$2j$
in figure 2
a). The structure of these waves is similar to the known stable elevation waves except for the local shape between two big troughs. However, the increase of the distance between the troughs weakens their coherence, therefore a ‘collision-like’ instability occurs. An example of this phenomenon is presented in figure 8, whereby the small perturbation accelerates the left trough and decelerates the right one. Figure 8(a) depicts the profiles of the perturbation, which is the most unstable normal mode with the corresponding growth rate
$\unicode[STIX]{x1D706}=0.00614$
. We present in (b,c) the comparisons between the numerical simulation of the full Euler equations and the prediction of the linear theory at
$t=600,750$
, which show great agreement. After
$t=750$
, the nonlinear effects become significant and a quasi-elastic collision is observed where both waves survive and the primary effect of the collision is a rapid phase shift (see the wave trajectories in 8
d and typical wave profiles in 8
e). The computation shown in figure 8 was carried out with
$L=300$
,
$N=8000$
and
$\unicode[STIX]{x0394}t=4\times 10^{-4}$
.
Figure 8. Instability of an elevation wave which consists of two depression waves connected by two small dimples. The steady solution, which corresponds to the point
$2A$
in figure 2(a) (
$\unicode[STIX]{x1D702}(0)=0.0425$
and
$c=1.39$
), was disturbed by its most unstable mode (
$\unicode[STIX]{x1D706}=0.00614$
) shown in (a) for
$\widetilde{\unicode[STIX]{x1D702}}$
(top) and
$\widetilde{\unicode[STIX]{x1D719}}$
(bottom). (b,c) Comparison between the linear theory (dashed line) and the nonlinear evolution code (solid line) at
$t=600$
and
$t=750$
respectively. (d) Time-dependent solution in the
$\unicode[STIX]{x1D709}-t$
plane. (e) Four free-surface profiles are shown from top to bottom at times
$t=750,1000,2000,8000$
(solid line), compared with the unperturbed steady solution (dotted line). The solution is shown in a frame of reference moving to the right with the speed of the undisturbed wave.
Lastly, we consider the stability of asymmetric gravity–capillary solitary waves. Following Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), the initial guess for finding an asymmetric solitary wave was obtained by superposing a depression wave (
$1a$
in figure 1
a) and a stable elevation one (
$2d$
in figure 2
a), both of which translate at
$c=1.4$
. With this initial guess, the Newton’s algorithm converged to the solution 9f in figure 9. It is natural to ask whether the obtained asymmetric solution is stable since its constituent parts are both stable.
Using 9f as the starting point, the bifurcation curve can be followed, which is shown in figure 9. The spectral stability analysis was then carried out using the method described above, and the growth rate is presented in table 2. We conclude from table 2 that all asymmetric waves on this branch are unstable, even when the wave profiles resemble two stable solitary waves placed side by side. It is noted that the investigation of the global bifurcation of asymmetric gravity–capillary solitary waves is beyond the scope of this paper, and the interested reader is referred to Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014).
Table 2. Instability growth rate
$\unicode[STIX]{x1D706}$
of asymmetric solitary waves for various values of the speed
$c$
. The energy of the wave is listed for reference.
Figure 9. Speed–energy bifurcation diagram of a branch of asymmetric solitary waves, together with typical wave profiles for
$c=1.4$
(9c and 9f) and
$c=1.39$
(9a and 9h). The branch arises from an asymmetric wave with the initial guess being composed of a stable elevation wave and a depression one propagating at
$c=1.4$
.
Figure 10. Time evolution of the unstable mode for an asymmetric solitary wave (9c in figure 9). All profiles are shown in a frame of reference moving to the right with
$c=1.4$
. (a) The profiles of the most unstable normal mode. (b,c) Comparison between the linear theory (dashed line) and the nonlinear computation (solid line) at
$t=125$
and
$t=160$
, respectively. (d) The evolution of the perturbation using the nonlinear evolution code. (e) Wave profiles at
$t=160,200,500,1200$
from top to bottom, with dotted line corresponding to the unperturbed steady solution.
Figures 10 and 11 show two numerical experiments of asymmetric solitary waves (9c and 9f respectively) perturbed with their most unstable normal modes, which were performed with
$L=300$
,
$N=8000$
and
$\unicode[STIX]{x0394}t=4\times 10^{-4}$
. Figure 10(b–d) show the nonlinear evolution of the initial disturbance whose profiles are given in figure 10(a). Figure 10(c) strongly suggests that the linear theory captures the development of the unstable mode until
$t\approx 160$
, a value of order
$1/\unicode[STIX]{x1D706}$
(here
$\unicode[STIX]{x1D706}=0.02608$
). Snapshots in figure 10(e) demonstrate the long-time behaviour of the instability, which manifests itself by a visible change of the central troughs and the subsequent breakdown of the coherent structure. Figure 11 depicts the instability of the wave 9f. The decrease of the amplitude of the elevation component increases its speed, while the increase of the amplitude of the depression component reduces its speed. Hence the small difference in speed between two components leads to changes of the connecting portion during the interaction. The comparisons between the nonlinear simulation based on (2.20) and the linear theory are made in figure 11(b,c), which show remarkable agreement. In contrast to the instability shown in figure 8, the elevation component is totally destroyed in the long-time dynamics.
Figure 11. Time-dependent computation of an asymmetric solitary waves (9f in figure 9) subject to the unstable normal mode. All profiles are shown in a frame of reference moving to the right with
$c=1.4$
. (a) Profiles of the most unstable normal mode. (b,c) Comparison between the linear theory (dashed line) and the nonlinear computation (solid line) at
$t=800$
and
$t=1000$
respectively. (d) From top to bottom, the snapshots show the evolution of the free surface at time
$t=1000,1500,2750,5000$
(solid line), and the unperturbed profile for reference.
4 Concluding remarks
The stability and dynamics of two-dimensional gravity–capillary solitary waves in deep water have been revisited numerically. The shedding of depression solitary waves and stable elevation waves have been achieved by moving a proper fully localised pressure distribution with the speed close to the phase speed minimum to mimic the jet of air impinging on the surface of a steady stream. The stability of depression solitary waves with an overhanging structure has been confirmed by Saffman’s theorem. While the numerical results that the overhanging structure survived the solitary-wave collisions provided further evidence of its stability.
The stability characteristics of multipacket solitary waves, including both symmetric and asymmetric waves, have been investigated using linear stability analysis and direct numerical simulations. Based on a simple rescaling, we have integrated hodograph transformation and time-dependent conformal map, so that the linear eigenmode obtained by the former technique can be used to prepare the initial data for the latter one. Comparisons between the time evolution of the unstable normal mode with the numerical integration of the full Euler equations and the prediction of the linear theory show excellent agreement for all solitary waves we have tested, validating the accuracy and reliability of our numerical scheme.
For symmetric elevation solitary waves, we have found that the stability exchange occurs twice along the bifurcation curve, which extends the previous study by Calvo & Akylas (Reference Calvo and Akylas2002) who were concerned only with the stability exchange near the first turning point. A special branch of asymmetric solitary waves has also been considered and found to be unstable, even when the wave profiles feature a depression wave and a stable elevation one separated by small oscillations. We have further found that some multipacket solitary waves exhibit collision-like behaviours in the long-time dynamics. Particularly, there is an oscillation before the separation of two solitary waves (see figure 8 d). This phenomenon is similar to the weak interaction of solitary waves in nonlinear optics where the governing equation is the generalised NLS equation (see, for example, figure 3 in Zhu & Yang (Reference Zhu and Yang2007) for cubic–quintic NLS equation). While beyond the scope of this paper, it is of great interest to study the similarities between hydrodynamics and nonlinear fibre optics.
Instabilities for some multipacket gravity–capillary solitary waves are weak (i.e. the growth rate is small), therefore they are probably observable in experiments within a reasonable timescale (1000 dimensionless time units correspond to approximately 17 s). It must be pointed out that due to the small typical wavelength (
${\approx}1.7~\text{cm}$
) of gravity–capillary waves, the viscous effect, which we have neglected, plays an important role even in the deep-water case. The Reynolds number, which measures the dissipation rate, is defined as
$$\begin{eqnarray}\displaystyle \mathit{Re}=\frac{\text{inertial forces}}{\text{viscous forces}}=\frac{\text{velocity scale}\times \text{length scale}}{\text{kinematic viscosity}}=\frac{1}{\unicode[STIX]{x1D708}}\left(\frac{\unicode[STIX]{x1D70E}^{3}}{\unicode[STIX]{x1D70C}^{3}g}\right)^{1/4}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D708}$
is the kinematic viscosity. In CGS units,
$g=981~\text{cm}~\text{s}^{-2}$
,
$\unicode[STIX]{x1D70E}=72.8~\text{g}~\text{s}^{-2}$
,
$\unicode[STIX]{x1D70C}=0.998~\text{g}~\text{cm}^{-3}$
and
$\unicode[STIX]{x1D708}=0.013~\text{cm}^{2}~\text{s}^{-1}$
for deep water. Hence the Reynolds number is approximately 340 for the air–water interface. The linear theory (Lamb (Reference Lamb1932), §§ 348–349) gives the rate of decay of order
$2/Re\approx 0.006$
, which is comparable to the growth rate of the most unstable modes of some solitary waves. To reduce the relative effect of viscosity in practice, using other fluids like mercury (see Falcon, Laroche & Fauve Reference Falcon, Laroche and Fauve2007 for example) or conducting the experiments in a microgravity environment is plausible.
Finally we make remarks about the stability of two-dimensional gravity–capillary waves to three-dimensional perturbations. Based on an asymptotic analysis, Kim & Akylas (Reference Kim and Akylas2007) showed that plane depression solitary waves in deep water are eventually unstable to transverse perturbations of sufficiently long wavelength. This fact is related to the existence of genuinely three-dimensional solitary waves, which are commonly referred to as ‘lumps’ in the literature. Gravity–capillary lumps were numerically found by Părău, Vanden-Broeck & Cooker (Reference Părău, Vanden-Broeck and Cooker2005) in the full Euler equations, and experimentally observed by Diorio et al. (Reference Diorio, Cho, Duncan and Akylas2009, Reference Diorio, Cho, Duncan and Akylas2011). The stability and dynamics of lumps, such as lump generation either by a localised forcing moving steadily or via transverse instability of a plane wave, long-time evolutions of unstable lumps, and collisions between stable lumps, have not been truly justified with the fully nonlinear water-wave equations in three dimensions, and are therefore of great interest for future studies.
Acknowledgements
The author would like to thank B. Akers (Air Force Institute of Technology) for asking the question that led to this work, and P. Milewski (University of Bath) for insightful discussions and suggestions on an early version of the paper. This work was supported by the National Natural Science Foundation of China (no. 11232012), and Key Research Program of Frontier Sciences, CAS (no. QYZDB-SSW-SYS015).
