1. Introduction
Hydroelasticity, a name adapted from aeroelasticity, is concerned with the motion and distortion of deformable bodies responding to hydrodynamic excitations, and the associated reactions on the motion of the environmental fluid. Hydroelastic waves enjoy wide usage in marine structures and sea transport. Modern applications of hydroelastic waves abound: very large floating structures usable as fully functional airport runways (Megafloat project in Japan); large fast merchant ships and container vessels that are relatively more flexible (Wang Reference Wang2000); flexible risers to transport hydrocarbon (mainly referring to oil) from the seabed to shore or offshore facilities (Jain Reference Jain1994); and safe use of lake and ocean ice for roadways and landing strips (Wilson Reference Wilson1958; Takizawa Reference Takizawa1985; Squire et al. Reference Squire, Robinson, Langhorne and Haskell1988). Due to this physical and industrial significance, in-depth knowledge of the characteristics of hydroelastic waves is therefore important.
The present work considers the irrotational motion of a two-dimensional inviscid and incompressible fluid of infinite depth, with the top surface being in contact with a frictionless thin elastic sheet. Since there are two restoring forces across the free surface (the gravity and the flexural rigidity due to the elastic bending) the hydroelastic waves propagating through the elastic cover are also called flexural–gravity waves. This problem has been proposed as an ideal model for the dynamics of a large body of waves covered by a floating ice sheet in polar regions, and has been studied extensively mostly in the linear setting since the pioneering work by Greenhill (Reference Greenhill1886) (the reader is referred to the monograph by Squire et al. (Reference Squire, Hosking, Kerr and Langhorne1996) which describes in detail the research based on the linear theory prior to 1996).
The linear theory is a good approximation when the deformations of the elastic cover are relatively small, but it becomes unreliable as the amplitude of the waves grows (the reader is referred to the reports of intense-in-ice events by Marko (Reference Marko2003) which highlight the limitations of the linear theory). In particular, the linear theory fails for the problem of a concentrated line load moving steadily over an ice sheet floating on a fluid when the velocity of the load is at the minimum of the phase speed (denoted by
$c_{min}$
). This breakdown of the linear theory is due to an accumulation of the energy which makes the displacement of the free surface grow without limit. In order to resolve the problem, Părău & Dias (Reference Părău and Dias2002) developed a nonlinear theory valid when the speed of the load is close to the critical value
$c_{min}$
. They conducted a normal form analysis which leads to a forced nonlinear Schrödinger equation (NLS) with the mean depth of the fluid being a parameter. Their results show that for water of sufficiently large depth, bounded responses in the form of hydroelastic solitary waves can exist for speeds up to
$c_{min}$
, while there is a range of forcing speed below
$c_{min}$
for which there are no steady solutions when the fluid is relatively shallow. Unsteady simulations using truncated models and high-order spectral methods were carried out by Bonnefoy, Meylan & Ferrant (Reference Bonnefoy, Meylan and Ferrant2009), confirming the predictions of Părău & Dias (Reference Părău and Dias2002) in the appropriate regime. Milewski, Vanden-Broeck & Wang (Reference Milewski, Vanden-Broeck and Wang2011) revisited the same problem and found hydroelastic solitary waves in deep water in the absence of moving loads, even though the weakly nonlinear analysis shows that the associated cubic NLS is of defocusing type at the minimum of the phase speed. However, these solitary waves are unusual, since they can only exist with finite amplitudes. In other words, they do not bifurcate from infinitesimal periodic waves like the gravity–capillary solitary waves in deep water (see, e.g., Vanden-Broeck & Dias Reference Vanden-Broeck and Dias1992; Wang, Vanden-Broeck & Milewski Reference Wang, Vanden-Broeck and Milewski2014) or from infinitesimal long waves like the Korteweg–de Vries (KdV) equation. Furthermore, the unsteady simulations presented in the same paper reveal that these solitary waves may arise naturally from the moving load problem subject to a moderate-amplitude near-critical forcing. This fact indicates that the existence and the stability of free solitary waves are crucial for an understanding of the forced problem.
All of the results mentioned in the last paragraph are based on the Kirchhoff–Love model, which uses
$\partial _{xx}{\it\kappa}$
to model the pressure due to the bending of the elastic sheet, where
${\it\kappa}$
is the curvature of the surface and
$x$
is the coordinate in the direction of wave propagation. The Kirchhoff–Love elastic model has also been used in the absence of moving loads, to study periodic hydroelastic waves (Forbes Reference Forbes1986, Reference Forbes1988), generalised flexural–gravity solitary waves (Vanden-Broeck & Părău Reference Vanden-Broeck and Părău2011) and the unsteady interaction between a fluid-loaded elastic plate and a mean flow (Peake Reference Peake2001). Although the Kirchhoff–Love model is widely used in the literature, it has some limitations. In particular, it does not appear to have an elastic potential. More recently, Toland (Reference Toland2007, Reference Toland2008) proposed a novel nonlinear elastic model using the Cosserat theory of hyperelastic shells satisfying Kirchhoff’s hypotheses, which has a clear variational structure. From then on, analytical and numerical investigations of this new model have been gradually carried out. Of note are the works of Toland (Reference Toland2008), who rigorously proved the existence of periodic hydroelastic waves, Guyenne & Părău (Reference Guyenne and Părău2012), who discovered that both elevation and depression branches exist below the minimum of the phase speed at finite amplitude in deep water, Wang, Vanden-Broeck & Milewski (Reference Wang, Vanden-Broeck and Milewski2013), who extended the branch of elevation solitary waves to the highly nonlinear regime with the wave profiles featuring multi-packet structure and computed periodic waves with an overhanging structure, Gao & Vanden-Broeck (Reference Gao and Vanden-Broeck2014), who investigated generalised solitary waves extensively, and Page & Părău (Reference Page and Părău2014), who considered nonlinear hydroelastic hydraulic falls past a submerged bottom obstruction.
There are relatively fewer studies on unsteady hydroelastic waves based on Toland’s model. It is worth mentioning that Guyenne & Părău performed direct numerical simulations for unsteady hydroelastic solitary waves in deep water (Guyenne & Părău Reference Guyenne and Părău2012) and shallow water (Guyenne & Părău Reference Guyenne and Părău2014). However, their numerics is based on reduced models with the truncated Dirichlet–Neumann operator, which cannot be used to study highly nonlinear waves, such as overhanging waves. To the best of our knowledge, there have been no computations of the dynamics of hydroelastic solitary waves using the full Euler equations and Toland’s model.
In this paper, we use a time-dependent conformal map technique to study numerically the bifurcation, stability, collision and excitation of hydroelastic solitary waves. We restrict our attention to deep water and to two dimensions. Motivated by the recent work of Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), we compute a new branch of solitary waves which are non-symmetric in the direction of wave propagation. These waves have a multi-packet structure. We study numerically the stability of all hydroelastic solitary waves found to date. Our results show that some symmetric solitary waves are robust subject to longitudinal perturbations, including not only moderate-amplitude single-trough depression solitary waves, but also certain elevation waves resembling two large troughs placed side by side. We show that these stable two-trough coherent structures can be excited by multiple loads simultaneously moving with a speed slightly below
$c_{min}$
.
The rest of the paper is structured as follows. We state the mathematical formulation and introduce the time-dependent conformal map technique used to reduce the dimension of the problem in § 2 and § 3 respectively. We review the properties of symmetric solitary waves in § 4.1, and compute asymmetric solitary waves in § 4.2. The stability properties for both symmetric and asymmetric waves are studied numerically in § 5. We then consider the interaction between two stable solitary waves, including both head-on and overtaking collisions. Numerical experiments on the excitation of solitary waves by one or multiple constant-velocity loads are performed in § 6.2. Concluding remarks and possible future projects are presented in the last section.
2. Formulation
We consider a two-dimensional irrotational flow of an incompressible and inviscid fluid beneath a thin elastic sheet. The fluid is assumed to be of infinite depth. We introduce Cartesian coordinates with the
$y$
-axis directed vertically upwards and with
$y=0$
at the undisturbed level of the elastic sheet. The deformation of the elastic sheet is denoted by
$y={\it\eta}(x,t)$
, where
$t$
is the time. Since the flow is irrotational we can introduce a potential function
${\it\phi}$
, such that the velocity field reads
$({\it\phi}_{x},{\it\phi}_{y})$
. The problem reduces then to solving Laplace’s equation
$$\begin{eqnarray}\displaystyle {\it\phi}_{xx}+{\it\phi}_{yy}=0,\quad \text{for }-\infty <y<{\it\eta}(x,t). & & \displaystyle\end{eqnarray}$$
The main approximations made here are that the elastic sheet is thin, and that its inertia and its stretching (or the existence of a prestressed state) are neglected. Therefore, the only restoring forces are gravity and the flexural elasticity of the sheet. In the present paper, we use the fully nonlinear model based on the special Cosserat theory of hyperelastic shells introduced by Toland (Reference Toland2008). Therefore, the nonlinear boundary conditions at
$y={\it\eta}(x,t)$
are the kinematic condition
$$\begin{eqnarray}\displaystyle {\it\eta}_{t}={\it\phi}_{y}-{\it\eta}_{x}{\it\phi}_{x} & & \displaystyle\end{eqnarray}$$
and the dynamic condition
$$\begin{eqnarray}\displaystyle {\it\phi}_{t}=-\frac{1}{2}\left[{\it\phi}_{x}^{2}+{\it\phi}_{y}^{2}\right]-g{\it\eta}-\frac{\mathscr{D}}{{\it\rho}}\mathscr{P}_{b}+\mathscr{P}_{e}(x,t), & & \displaystyle\end{eqnarray}$$
where
$g$
is the acceleration due to gravity,
${\it\rho}$
is the density of the fluid,
$\mathscr{P}_{e}$
is the external pressure distribution exerted on the elastic sheet,
$\mathscr{P}_{b}$
is the pressure distribution due to elastic bending and
$\mathscr{D}$
is the coefficient of flexural rigidity, defined as
$\mathscr{D}=Eh^{3}/12(1-{\it\nu}^{2})$
. Here,
$E$
is the Young’s modulus,
${\it\nu}$
is the Poisson ratio and
$h$
is the thickness of the elastic sheet. It follows from Toland (Reference Toland2008) (equation (1.13) in Toland’s paper when the surface tension is neglected) that
$$\begin{eqnarray}\displaystyle \mathscr{P}_{b}={\it\kappa}_{ss}+{\textstyle \frac{1}{2}}{\it\kappa}^{3}, & & \displaystyle\end{eqnarray}$$
where
${\it\kappa}$
is the curvature of the free surface and
$s$
is the arclength parameter. If the wave profile is single-valued, the curvature can be expressed in Cartesian coordinates as
$$\begin{eqnarray}\displaystyle {\it\kappa}=\frac{{\it\eta}_{xx}}{(1+{\it\eta}_{x}^{2})^{3/2}}. & & \displaystyle\end{eqnarray}$$
Finally, the far-field condition
$$\begin{eqnarray}\displaystyle {\it\phi}_{y}\rightarrow 0\quad \text{as }y\rightarrow -\infty & & \displaystyle\end{eqnarray}$$
completes the formulation. The system can be non-dimensionalised by choosing
$$\begin{eqnarray}\displaystyle \left[\frac{\mathscr{D}}{{\it\rho}g}\right]^{1/4},\quad \left[\frac{\mathscr{D}}{{\it\rho}g^{5}}\right]^{1/8},\quad \left[\frac{g\mathscr{D}^{3}}{{\it\rho}^{3}}\right]^{1/8} & & \displaystyle\end{eqnarray}$$
as the units of length, time and potential respectively. The dynamic boundary condition can be rewritten as
$$\begin{eqnarray}\displaystyle {\it\phi}_{t}=-{\textstyle \frac{1}{2}}[{\it\phi}_{x}^{2}+{\it\phi}_{y}^{2}]-{\it\eta}-\left({\textstyle \frac{1}{2}}{\it\kappa}^{2}+{\it\kappa}_{ss}\right)+\mathscr{P}_{e}(x,t). & & \displaystyle\end{eqnarray}$$
Equations (2.1), (2.2), (2.5), (2.6) and (2.8) form a Hamiltonian system, with the action functional being the total energy of the fluid, which is the sum of the kinetic energy and the potential energy:
$$\begin{eqnarray}\displaystyle \mathscr{E}=\frac{1}{2}\left[\int _{\mathbb{R}}\,\text{d}x\int _{-\infty }^{{\it\eta}}({\it\phi}_{x}^{2}+{\it\phi}_{y}^{2})\,\text{d}y+\int _{\mathbb{ R}}{\it\eta}^{2}\,\text{d}x+\int \,{\it\kappa}^{2}\,\text{d}s\right]. & & \displaystyle\end{eqnarray}$$
Here,
$\mathbb{R}$
is the real line. We denote the velocity potential on the free surface by
${\it\varphi}(x,t)\triangleq {\it\phi}(x,{\it\eta}(x,t),t)$
. Working in the canonical variables
${\it\varphi}$
and
${\it\eta}$
, the kinematic and dynamic boundary conditions can be recast as
$$\begin{eqnarray}\displaystyle {\it\eta}_{t}=\frac{{\it\delta}\mathscr{E}}{{\it\delta}{\it\varphi}},\quad {\it\varphi}_{t}=-\frac{{\it\delta}\mathscr{E}}{{\it\delta}{\it\eta}}. & & \displaystyle\end{eqnarray}$$
It is noted that the effects due to the inertia and the thickness of the elastic sheet can also be incorporated into the pressure equation (2.8) (see Squire et al. (Reference Squire, Hosking, Kerr and Langhorne1996) for the expression of the inertia and Forbes (Reference Forbes1986) for how to model the system with the thickness of the plate). However, the unsteady simulations become much more complicated when the inertia is taken into account. For simplicity, we neglect these two effects, and only consider the pressure jump exerted by the elastic sheet due to flexing. Even with these approximations, the hydroelastic wave problem has been modelled with a variety of methods. The reader is referred to Milewski & Wang (Reference Milewski and Wang2013) for a discussion.
3. Numerical method
The main idea for solving two-dimensional full Euler equations is based on a time-dependent conformal mapping. This numerical method was pioneered by Dyachenko et al. (Reference Dyachenko, Kuznetsov, Spector and Zakharov1996). We start the derivation with finding a transformation that maps the physical domain
$-\infty <y<{\it\eta}(x,t)$
into the lower half-plane with horizontal and vertical coordinates denoted by
${\it\xi}$
and
${\it\zeta}$
respectively. The conformal mapping from
$(x,y)$
to
$({\it\xi},{\it\zeta})$
can be found by solving the following boundary-value problem:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}ll@{}}\displaystyle y_{{\it\xi}{\it\xi}}+y_{{\it\zeta}{\it\zeta}}=0 & \text{for }-\infty <{\it\zeta}<0,\\ \displaystyle y=Y({\it\xi},t) & \text{at }{\it\zeta}=0,\\ \displaystyle y\sim {\it\zeta} & \text{as }{\it\zeta}\rightarrow -\infty ,\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$Y({\it\xi},t)={\it\eta}(x({\it\xi},0,t),t)$
. The harmonic conjugate
$x({\it\xi},{\it\zeta},t)$
is defined through the Cauchy–Riemann relations for the complex function
$x({\it\xi},{\it\zeta},t)+\text{i}y({\it\xi},{\it\zeta},t)$
. In the transformed plane, the velocity potential
${\it\phi}$
and its harmonic conjugate
${\it\psi}({\it\xi},{\it\zeta},t)$
also satisfy Laplace’s equations. Defining
${\it\Phi}({\it\xi},t)\triangleq {\it\phi}(x({\it\xi},0,t),y({\it\xi},0,t),t)$
,
${\it\Psi}({\it\xi},t)\triangleq {\it\psi}({\it\xi},0,t)$
and
$X({\it\xi},t)\triangleq x({\it\xi},0,t)$
, we obtain after some elementary analysis
$$\begin{eqnarray}\displaystyle X={\it\xi}-\mathscr{H}[Y]\quad \text{and}\quad {\it\Psi}=\mathscr{H}[{\it\Phi}], & & \displaystyle\end{eqnarray}$$
where
$\mathscr{H}$
is the operator of Hilbert transformation with the Fourier symbol
$\text{i}\text{sgn}(k)$
. It can also be defined in the physical space as
$$\begin{eqnarray}\displaystyle \mathscr{H}[f]({\it\xi})=\text{PV}\int _{\mathbb{R}}\frac{f({\it\xi}^{\prime })}{{\it\xi}^{\prime }-{\it\xi}}\,\text{d}{\it\xi}^{\prime }, & & \displaystyle\end{eqnarray}$$
where ‘PV’ indicates the Cauchy principal value of the integral. By using the chain rule, the kinematic boundary condition in the transformed plane takes the form
$$\begin{eqnarray}\displaystyle X_{{\it\xi}}Y_{t}-Y_{{\it\xi}}X_{t}=-{\it\Psi}_{{\it\xi}}. & & \displaystyle\end{eqnarray}$$
Following Dyachenko et al. (Reference Dyachenko, Kuznetsov, Spector and Zakharov1996), the evolution equation for
$Y$
is
$$\begin{eqnarray}\displaystyle Y_{t}=Y_{{\it\xi}}\mathscr{H}\left[\frac{{\it\Psi}_{{\it\xi}}}{J}\right]-X_{{\it\xi}}\frac{{\it\Psi}_{{\it\xi}}}{J}, & & \displaystyle\end{eqnarray}$$
where
$J=X_{{\it\xi}}^{2}+Y_{{\it\xi}}^{2}$
is the Jacobian of the conformal map. Finally, we reformulate the dynamic boundary condition using the new variable
${\it\xi}$
as
$$\begin{eqnarray}\displaystyle {\it\Phi}_{t}=\frac{{\it\Psi}_{{\it\xi}}^{2}-{\it\Phi}_{{\it\xi}}^{2}}{2J}-Y-\frac{1}{2}\left[\frac{{\it\kappa}_{{\it\xi}{\it\xi}}}{J}+\left(\frac{{\it\kappa}_{{\it\xi}}}{J}\right)_{{\it\xi}}+{\it\kappa}^{3}\right]+{\it\Phi}_{{\it\xi}}\mathscr{H}\left[\frac{{\it\Psi}_{{\it\xi}}}{J}\right]+\mathscr{P}_{e}, & & \displaystyle\end{eqnarray}$$
where the curvature is now written as
$$\begin{eqnarray}\displaystyle {\it\kappa}=\frac{X_{{\it\xi}}Y_{{\it\xi}{\it\xi}}-X_{{\it\xi}{\it\xi}}Y_{{\it\xi}}}{J^{3/2}}. & & \displaystyle\end{eqnarray}$$
In order to find the dispersion relation of the system, we linearise the surface Euler system (3.5) and (3.6) by taking
$Y$
,
${\it\Phi}_{{\it\xi}}$
,
${\it\Psi}_{{\it\xi}}$
small,
$\mathscr{P}_{e}=0$
and
$X_{{\it\xi}}\sim 1$
,
$J\sim 1$
. This yields
${\it\Phi}_{tt}=\mathscr{H}[{\it\Psi}_{{\it\xi}}]+\mathscr{H}[{\it\Psi}_{{\it\xi}{\it\xi}{\it\xi}{\it\xi}{\it\xi}}]$
. Therefore, the dispersion relation is
$$\begin{eqnarray}\displaystyle {\it\omega}^{2}=k^{2}c^{2}=|k|(1+k^{4}). & & \displaystyle\end{eqnarray}$$
Relation (3.8) implies that the phase velocity
$c$
reaches its minimum value
$c_{min}\approx 1.3247$
when
$k=(1/3)^{1/4}$
. The phase velocity is then equal to the group velocity at this minimum.
We seek fully localised solutions travelling with speed
$c$
. In order to guarantee that there are no waves in the far field, we choose
$c<c_{min}$
. We assume that the functions
$Y$
and
${\it\Phi}$
depend on
${\it\xi}-ct$
. It then follows from (3.4) that
$$\begin{eqnarray}\displaystyle {\it\Psi}_{{\it\xi}}=cY_{{\it\xi}},\quad \text{and hence}\quad {\it\Phi}_{{\it\xi}}=-c\mathscr{H}[Y_{{\it\xi}}]. & & \displaystyle\end{eqnarray}$$
From (3.5), one can conclude that
$$\begin{eqnarray}\displaystyle \mathscr{H}\left[\frac{{\it\Psi}_{{\it\xi}}}{J}\right]=c\left(\frac{X_{{\it\xi}}}{J}-1\right). & & \displaystyle\end{eqnarray}$$
Substituting (3.9) and (3.10) into the dynamic boundary condition (3.6) with
$\mathscr{P}_{e}=0$
(since we seek free solitary waves) and noticing that
${\it\Phi}_{t}=-c{\it\Phi}_{{\it\xi}}$
, we obtain after some algebra the following nonlinear equation:
$$\begin{eqnarray}\displaystyle \frac{c^{2}}{2}\left(\frac{1}{J}-1\right)+Y+\frac{1}{2}\left[\frac{{\it\kappa}_{{\it\xi}{\it\xi}}}{J}+\left(\frac{{\it\kappa}_{{\it\xi}}}{J}\right)_{{\it\xi}}+{\it\kappa}^{3}\right]=0. & & \displaystyle\end{eqnarray}$$
This, together with
$X_{{\it\xi}}=1-\mathscr{H}[Y_{{\it\xi}}]$
, defines an integro-differential system. Finally, the velocity potential can be recovered from
${\it\Phi}=-c\mathscr{H}[Y]$
. We shall emphasise that the curvature (2.5) is only valid for single-valued solutions, while the formulation (3.11) based on the new parameter
${\it\xi}$
is even true for multivalued wave profiles, and the system (3.5), (3.6) allows us to compute the evolution of overhanging waves. In the new coordinates, the Hamiltonian (2.9) becomes
$$\begin{eqnarray}\displaystyle \mathscr{E}=\frac{c^{2}}{2}\int _{\mathbb{R}}Y_{{\it\xi}}\mathscr{H}[Y]\,\text{d}{\it\xi}+\frac{1}{2}\int _{\mathbb{R}}Y^{2}X_{{\it\xi}}\,\text{d}{\it\xi}+\frac{1}{2}\int _{\mathbb{R}}\frac{(Y_{{\it\xi}{\it\xi}}X_{{\it\xi}}-X_{{\it\xi}{\it\xi}}Y_{{\it\xi}})^{2}}{J^{5/2}}\,\text{d}{\it\xi}, & & \displaystyle\end{eqnarray}$$
where the terms on the right-hand side correspond to kinetic energy, gravitational potential energy and elastic potential energy respectively. We remark that conformal mapping is a conventional technique to handle free-surface water wave problems in two dimensions. Equation (3.11) and its modified versions have been widely used in computing periodic and solitary waves in deep water. Of note are the works of Crapper (Reference Crapper1957), who obtained analytical solutions for pure capillary waves in terms of elementary functions, Schwarts & Vanden-Broeck (Reference Schwarts and Vanden-Broeck1979) and Longuet-Higgins (Reference Longuet-Higgins1988), who computed different branches of periodic capillary–gravity waves, and Guyenne & Părău (Reference Guyenne and Părău2012), who considered flexural–gravity solitary waves based on Toland’s nonlinear elastic model. All of the aforementioned works show that multivalued steady profiles exist in capillary/capillary–gravity/flexural–gravity waves, with a limiting configuration pitching off a closed bubble.
Details of the numerical procedure for finding fully localised travelling waves in (3.11) are given as follows. The solitary-wave solutions are approximated by long periodic waves. Therefore, (3.11) can be solved numerically via truncating the Fourier series of
$Y$
, namely
$$\begin{eqnarray}\displaystyle Y({\it\xi})=\mathop{\sum }_{n=0}^{N}a_{n}\cos (n{\rm\pi}{\it\xi}/L)+\mathop{\sum }_{n=1}^{N}\sin (n{\rm\pi}{\it\xi}/L). & & \displaystyle\end{eqnarray}$$
The Fourier coefficients
$a_{n}$
and
$b_{n}$
need to be found so that the dynamic boundary condition (3.11) is satisfied. This is done numerically. We first discretise the domain
$[-L,L)$
into a uniform mesh. We evaluate (3.11) on the grid points and then project it onto each element of the
$\cos (n{\rm\pi}{\it\xi}/L)$
and
$\sin (n{\rm\pi}{\it\xi}/L)$
basis for
$n=0,1,\ldots ,N$
. All of the derivatives and the Hilbert transform are calculated via Fourier multipliers, making the programme efficient and accurate, while the nonlinear terms are computed in real space. The resultant system of nonlinear algebraic equations is solved by Newton’s method. The underlying period
$2L$
and the total wavenumber
$N$
are both chosen to be sufficiently large so that the wave profiles hardly change as
$L$
and
$N$
are further increased. In most computations, we keep the mesh size
${\rm\Delta}{\it\xi}=0.05$
. We stop the iterations when the
$l^{\infty }$
-norm of the residual error is less than
$10^{-10}$
. It is noted that for symmetric solitary waves
$b_{n}=0$
for all
$n$
, which reduces the computing time considerably. As a validation of the numerical method, we recomputed the solutions of Crapper (Reference Crapper1957) for periodic capillary waves of infinite depth. We obtained an excellent agreement with Crapper’s solutions, and this was tested by checking that the
$l^{\infty }$
-norms of the errors were all less than
$10^{-14}$
. In particular, we were able to compute the overhanging profiles which occur for waves close to the limiting configuration.
4. Travelling waves of permanent form
4.1. Symmetric waves
In this section we describe the basic properties of symmetric solitary waves propagating at a constant velocity
$c$
without change of form. These waves will be used in § 4.2 as building blocks to construct new asymmetric waves. The computation of steady symmetric solitary waves can also serve as a validation of the numerical method since we just reproduce the results in Guyenne & Părău (Reference Guyenne and Părău2012) and Wang et al. (Reference Wang, Vanden-Broeck and Milewski2013).
Figure 1. Speed–amplitude and speed–energy curves for symmetric elevation and depression solitary waves emerging from the bifurcation point
$c^{\ast }\approx 1.3247$
. (a) Speed–amplitude curves of the elevation and depression branches. The value of
${\it\eta}$
at the middle point is considered as the amplitude. The depression branch is monotonic for
$c\in [0,c^{\ast })$
(only part of the curve is shown), while the elevation branch demonstrates a complex behaviour with multiple turning points. (b) Speed–energy curve of the elevation branch, showing a zigzag behaviour. (c) Speed–energy curve of the depression branch with two stationary points. The energy partition is shown in (d) for the elevation branch and in (e) for the depression branch, where the total energy is partitioned into kinetic energy (solid line), gravitational potential energy (dotted line) and elastic potential energy (dash-dotted line).
Guyenne & Părău (Reference Guyenne and Părău2012) first computed free solitary waves for the Toland model. They found that there is an elevation branch (a positive free-surface elevation
${\it\eta}(0)$
at the centre of the wave) and a depression branch (a negative free-surface elevation
${\it\eta}(0)$
at the centre of the wave) bifurcating from the phase speed minimum
$c^{\ast }\approx 1.3247$
. They noticed that these solitary waves exist only with non-zero amplitude, as shown in figure 1(a). A similar result was found earlier for the Kirchhoff–Love model (see Milewski et al.
Reference Milewski, Vanden-Broeck and Wang2011). Later, Wang et al. (Reference Wang, Vanden-Broeck and Milewski2013) calculated the complete solution branches for depression and elevation solitary waves. They found that the elevation branch exhibits multiple turning points, which indicates that different elevation solitary waves can travel at the same speed (a typical profile of elevation solitary waves is shown in figure 7
b). For the depression branch, both Guyenne & Părău (Reference Guyenne and Părău2012) and Wang et al. (Reference Wang, Vanden-Broeck and Milewski2013) found that as the speed
$c$
decreases and approaches zero, the wave profiles become steeper, eventually with an overhanging structure (a typical example of overhanging depression solitary waves is shown in figure 6
b).
The two branches in figure 1(a) bifurcate from the phase speed minimum
$c^{\ast }\approx 1.3247$
with non-zero amplitudes. This is to be contrasted with the problem of gravity–capillary waves, where the branches bifurcate from the phase speed minimum with
${\it\eta}(0)=0$
(see, for example, Vanden-Broeck & Dias Reference Vanden-Broeck and Dias1992). The reason for this difference is that the associated NLS is of defocusing type for the present problem (see Milewski & Wang Reference Milewski and Wang2013 for details).
The speed–energy bifurcation curves shown in figure 1(b,c) are closely related to the stability properties of solitary waves subject to longitudinal perturbations. Saffman (Reference Saffman1985) considered the stability of periodic waves due to superharmonic perturbations (the perturbations are periodic in one wavelength of the basic wave), and gave a necessary but not sufficient condition for stability exchange. He pointed out that stability exchanges can only occur at critical points, either stationary points or turning points, of the speed–energy bifurcation curve, namely
$$\begin{eqnarray}\frac{\partial \mathscr{E}}{\partial c}=0\quad \text{or}\quad \frac{\partial c}{\partial \mathscr{E}}=0.\end{eqnarray}$$
Saffman’s argument is based on the Hamiltonian formulation of water waves (see Zakharov Reference Zakharov1968) and is developed for pure gravity waves. However, the result can be generalised to hydroelastic waves without any essential modification due to the Hamiltonian structure (2.9). Furthermore, since all of the disturbances are superharmonic for solitary waves, Saffman’s conclusion is extremely useful for our problem. In § 5, we will focus on the stability properties for solitary waves via direct numerical simulations together with Saffman’s result. In figure 1(d,e), the total energy is partitioned into three parts: kinetic energy (solid line), gravitational potential energy (dotted line) and elastic potential energy (dash-dotted line). We note that for the elevation branch the different types of energy are similar in their behaviour but differ quantitatively, namely kinetic
$\text{energy}>\text{gravitational}$
potential
$\text{energy}>\text{elastic}$
potential energy (see figure 1
d). For large-amplitude depression solitary waves (figure 1
e), as the translating speed decreases, the kinetic energy decreases and vanishes at
$c=0$
, while the elastic potential energy increases and reaches its global maximum at the static state.
4.2. Asymmetric waves
Asymmetric gravity–capillary waves in deep water were recently computed by Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014) for the full Euler equations. Their existence is related to the possibility of having several symmetric waves travelling at the same speed. Figure 1(a) shows that for values of
$c$
close to 1.3247, we also have several symmetric flexural–gravity waves travelling at the same speed. This suggests that we search for asymmetric waves for the present problem. In this section we provide numerical evidence that such waves do exist.
Figure 2. An initial guess for the computation of an asymmetric hydroelastic solitary wave (c), which is composed of a depression wave (a) and a one-hump elevation wave (b), both propagating at
$c=1.32$
.
In order to compute asymmetric hydroelastic solitary waves, it is essential to choose a good initial guess for Newton’s iteration. Following the procedure described in Wang et al. (Reference Wang, Vanden-Broeck and Milewski2014), we choose a depression wave (see figure 2
a) and an elevation wave (see figure 2
b) computed in the last section. These two waves are chosen so that they travel at the same speed
$c=1.32$
(this choice is possible for values of
$c$
close to the phase speed minimum). We shift the profiles and then merge them in the middle as a new wave profile (see figure 2
c). This new profile is considered as the initial guess of Newton’s method, and the algorithm converges to a solution after several iterations. After obtaining a convergent solution, we can follow the branch by using a continuation method. This yields the speed–energy bifurcation diagram shown in figure 3.
Figure 3. Speed–energy bifurcation diagram of asymmetric hydroelastic solitary waves. The solid curves correspond to the branch of asymmetric waves. The dotted curve and the dashed curve correspond to two branches of symmetric waves where the asymmetric branch bifurcates from. The labels (a–f) correspond to those waves travelling at speed
$c=1.31$
, whose profiles are shown in figure 4.
We present in figure 3 values of the energy of asymmetric waves versus the speed
$c$
. The solid curves correspond to asymmetric waves, whereas the dotted curve and the dashed curve correspond to symmetric waves. Typical wave profiles corresponding to the points labelled (a–h) are shown in figure 4, including two asymmetric ones (figure 4
c,d). The curves of figure 3 show that the speed of asymmetric solitary waves is always below the minimum of the phase speed. Therefore, these waves do not resonate with linear periodic waves, namely they do not turn into generalised solitary waves or into periodic waves (see Milewski et al.
Reference Milewski, Vanden-Broeck and Wang2011 for a discussion). As can be seen from figure 3, the solid branches join the dashed and the dotted branches, i.e. the asymmetric waves eventually become symmetric (see profiles (g,h) in figure 4). This shows that the asymmetric waves appear from a spontaneous symmetry-breaking bifurcation and vanish at another symmetry-breaking bifurcation. It is worth noting that profiles 4(a,b) from the upper branches essentially consist of two identical elevation waves, while the profiles 4(e,f) from the lower branches are two depression waves merged at the origin. Energywise, the results are quite reasonable since an elevation wave possesses more energy than a depression wave with the same speed.
Figure 4. Typical profiles corresponding to points (a–h) indicated in figure 3: (c,d) are the typical profiles of asymmetric waves; (g,h) correspond to the symmetry-breaking bifurcation points. All profiles are shown in the physical space.
There exist many other families of asymmetric solitary-wave solutions which can be found by using different initial guesses because of the many possible choices of elevation waves travelling at the same speed, as demonstrated in figure 1. We leave this as a future research interest and focus on the stability problem.
5. Stability of solitary waves
In this section, we investigate the stability of the waves computed in § 4 when they are subjected to longitudinal perturbations. This is achieved by solving numerically the evolution equations (3.5) and (3.6) with a (slightly perturbed) travelling wave as the initial condition. Both symmetric and asymmetric waves are considered. For convenience we choose a frame of reference moving with the speed of the unperturbed solitary wave.
5.1. Symmetric waves
We start with symmetric depression solitary waves. The speed–energy bifurcation curve for the depression branch shown again in figure 5 has two stationary points labelled (1) and (2). According to Saffman (Reference Saffman1985), stability exchanges may occur at these points. We first investigate the stability of the waves before the first stationary point (
$c\lesssim 0.8$
). The waves in this regime are large and steep, and some are even overhanging. Therefore, the numerics becomes very stiff if an explicit time integration method is used because of the strong constraint on the time step. We remove this constraint by using the backward Euler’s method which is implicit. We take an overhanging depression solitary wave
${\it\eta}({\it\xi})$
with
$c=0.5$
and
${\it\eta}(0)=-2.70$
. At
$t=0$
, a small perturbation
$0.01\cos ({\it\xi}){\it\eta}({\it\xi})$
is added to the steady solution. The system (3.5) and (3.6) is then integrated by the backward Euler’s method with the time step
$\text{d}t=0.002$
, and at each step the discretised nonlinear algebraic system is solved by Newton’s method. The time evolution of the perturbed wave is shown in figure 6(a). The profiles in the physical space at time
$t=0,2,4$
are presented in figure 6(b–d). We can conclude that the bubble changes its shape quickly and then becomes unstable. We stop the computation when the profile touches itself, forming a closed bubble. Other depression solitary waves lying in the segment
$0<c\lesssim 0.8$
are all found to be unstable based on similar numerical experiments.
Figure 5. Stability of the depression branch. Waves from the solid segment are stable whereas those from the dash-dotted parts are unstable. The stationary points have been marked as asterisks on the graph.
Figure 6. Time evolution of a large-amplitude depression solitary wave propagating with the velocity
$c=0.5$
and amplitude
${\it\eta}(0)=-2.70$
. The computation was performed with
$\text{d}t=0.002$
and
$\text{d}{\it\xi}=0.01$
. (a) Time evolution of the perturbation which is initially
$0.01\cos ({\it\xi}){\it\eta}({\it\xi})$
. (b–d) Wave profiles in the physical space at time
$t=0,2,4$
. Only the main parts of the solutions are shown in the profiles.
For the remaining time-evolution computations, we use a fourth-order Runge–Kutta method rather than the backward Euler’s scheme since the problems are not extremely stiff for moderate-amplitude waves. In most of the computations, the mesh size in space and the step size in time are chosen as
$\text{d}{\it\xi}=0.1$
and
$\text{d}t=5\times 10^{-4}$
respectively.
For the solid part in figure 5 (the segment between the stationary points (1) and (2)), a variety of perturbations with 5 % of the energy of the initial depression solitary waves do not show instability. However, as the speed increases and passes the second stationary point (2), the solitary wave becomes unstable again. Snapshots of the dynamics of a depression solitary wave in this region subject to a 1 % amplitude-decreasing perturbation are shown in figure 7(a,c,e). The wave eventually disperses out as time evolves. Finally, we can draw the conclusion that for the depression branch of hydroelastic solitary waves, there are two stationary points in the speed–energy bifurcation diagram, and that stability exchanges occur at both of them.
Figure 7. Time evolution of a perturbed depression wave (
$c=1.32$
) at
$t=0$
, 500, 1000 (a,c,e) and a perturbed elevation wave (
$c=1.323$
) at
$t=0$
, 1500, 3000 (b,d,f).
The stability problem for the elevation branch is more complicated since the speed–energy curve has many stationary points and turning points, and we need to choose one typical example between every two successive critical points and check its stability. In the speed–energy bifurcation diagram (figure 8), the stationary points and turning points have been numbered from
$\unicode[STIX]{x2460}$
to
$\unicode[STIX]{x2466}$
. We found that only the waves between the turning point
$\unicode[STIX]{x2461}$
and the stationary point
$\unicode[STIX]{x2462}$
(solid curve in figure 8) are stable when subjected to longitudinal perturbations. All other waves on the elevation branch (dot-dashed curve in figure 8) turn out to be unstable.
Figure 8. Speed–energy bifurcation curve for the branch of elevation; (b) is a zoom in of the box in (a). Waves from the solid branch are stable whereas those from the dot-dashed curves are unstable. The stationary points and the turning points are marked as asterisks and pentagrams respectively.
We present some time evolutions of elevation waves to demonstrate their stability properties. A wave from the right of point
$\unicode[STIX]{x2460}$
has been perturbed by 1 % of the amplitude of the steady wave. The snapshots show a small focusing phenomenon and a symmetry-breaking instability as well, and finally become a stable depression solitary wave with a radiated wave field of linear wavepackets shedding on both sides (see figure 7
b,d,f). Figure 9(a,c,e,g) displays a stable elevation wave which is located between points
$\unicode[STIX]{x2461}$
and
$\unicode[STIX]{x2462}$
of figure 8. These stable elevation solitary waves feature two large troughs separated by a small dimple. A 5 % amplitude-decreasing perturbation has been applied to the initial solitary wave but it does not show any instability for a long-time computation. It is worth mentioning that the amplitude-decreasing perturbation results in a wave of slightly smaller energy but larger speed than the original one, which causes the right translation in the moving frame. However, not all elevation waves featuring two depression solitary waves placed side by side are stable. Figure 9(b,d,f,h) presents an example of unstable elevation waves composed of two troughs. The wave is located on the very last branch, i.e. beyond
$\unicode[STIX]{x2466}$
, with
$c=1.25$
,
${\it\eta}(0)=0.113$
and
$\mathscr{E}=9.34$
. We choose a specific perturbation which decreases the amplitude of the left trough by 0.5 % and increases the right one by 0.5 %. Therefore, the speed of the left trough is greater than that of the right trough so that an overtaking collision occurs, from which we can deduce the instability.
Figure 9. Time evolution of a perturbed elevation wave (
$c=1.2$
) at
$t=0$
, 3000, 4000, 5000 (a,c,e,g) and a perturbed elevation wave (
$c=1.25$
) at
$t=0$
, 750, 1000, 1500 (b,d,f,h).
5.2. Asymmetric waves
We now move on to the stability of asymmetric solitary waves. It turns out that the whole branch computed in § 4.2 is unstable. One typical example is shown in figure 10, where we take an asymmetric solitary wave featuring a depression wave and an elevation wave connected with small ripples as the initial data (figure 10 a). As time goes on, the trough further to the left remains unchanged but the elevation part on the right abruptly changes to a large trough (figure 10 b). The right trough travels leftwards since its speed is less than the speed of the moving frame (figure 10 c). Therefore, the stability problem becomes an overtaking collision and only one large trough survives (figure 10 d). Numerical experiments have been performed for other asymmetric solitary waves on the same branch. These were found to be unstable in all cases tried.
Figure 10. Time evolution of an asymmetric solitary wave (
$c=1.3$
) perturbed by numerical noise. Snapshots are shown for
$t=0$
, 250, 750, 1250 from (a) to (d) respectively.
6. Dynamics of solitary waves
In this section we supplement the stability analysis of § 5 by studying the interaction of solitary waves and their generation by moving disturbances.
6.1. Collision
We first consider head-on collisions and overtaking collisions between stable hydroelastic solitary waves. In all of the collision experiments, the initial data were constructed by first shifting the existing travelling waves so as to minimise their overlap, and then adding them together.
A typical example of a head-on collision is shown in figure 11. Here, a stable elevation wave travelling rightwards collides with a stable depression wave travelling leftwards. The collision is inelastic even though both waves survive after the collision, since many visible ripples are shed on both sides during the interaction. Several other computations of head-on collisions were carried out, and in all cases the collisions appeared to be inelastic.
Figure 11. Head-on collision. An elevation wave travels rightwards at
$c=1.3$
and a depression wave travels leftwards at
$c=1.3$
. The wave profiles are shown for time
$t=0$
, 15, 30, 45, 60, 75, 90.
The second set of numerical experiments on the dynamics of hydroelastic solitary waves is on overtaking collisions. We choose two stable depression waves travelling in the same direction but with different translating speeds. When the amplitude difference between the two is large, only one wave survives after the collision, with a radiation field. However, both depression waves can survive if their amplitude difference is small. As can be seen from figure 12(a), the overtaking collision is also inelastic when both waves survive, since some radiation field is generated during the interaction. The accuracy of the integration was tested by monitoring the value of the total energy, the error of which was of order
$10^{-12}$
through the whole computation. The difference between the initial data and the solution obtained by reversing time after
$t=4200$
in the collision was of order
$10^{-9}$
(see figure 12
b). This further validates the time-dependent codes. We can also choose a stable elevation wave to overtake a depression one, and vice versa. One example of this kind of collision is presented in figure 13. The structure of the elevation solitary wave with two well-separated troughs was destroyed during the interaction, and two depression waves travelling at different speeds were ultimately observed.
Figure 12. Overtaking collision. Two depression waves travel rightforwards at
$c=1.26$
and
$c=1.25$
respectively. (a) Snapshots for time
$t=0$
, 700, 1400, 2100, 2800, 3500, 4200 are shown in a frame of reference moving with the speed
$c=1.25$
. (b) The difference between the initial data and the time-reversed solution computed from
$t=4200$
back to
$t=0$
.
Figure 13. Overtaking collision. An elevation wave travels rightforwards at
$c=1.3$
and a depression travels at
$c=1.25$
in the same direction. Snapshots for time
$t=0$
, 400, 800, 1200, 1600, 2000 are shown in a frame of reference moving with the speed
$c=1.25$
.
Figure 14. Snapshots of the free surface due to large forcing moving pressures (
$A=0.2$
and
$u=1.3$
) at
$t=25$
, 100, 140, 200, 325 (a–e). The force is switched on at
$t=0$
and off at
$t=125$
. It is initially placed at
$x=-220$
.
6.2. Excitation
The coupling between a moving load and the hydroelastic waves generated has attracted interest for a long time because of its applications in polar regions (Wilson Reference Wilson1958; Takizawa Reference Takizawa1985; Squire et al. Reference Squire, Robinson, Langhorne and Haskell1988). According to the linear theory, the critical velocity for resonance between the moving load and the hydroelastic waves occurs at the minimum of the phase speed. The resonance results in an accumulation of energy and an unlimited growth of the elastic-sheet displacement. Therefore, the linear theory fails to predict the physical behaviour of the system, and a nonlinear theory needs to be developed. This was considered by Părău & Dias (Reference Părău and Dias2002), Bonnefoy et al. (Reference Bonnefoy, Meylan and Ferrant2009) and Milewski et al. (Reference Milewski, Vanden-Broeck and Wang2011) for the Kirchhoff–Love model and by Guyenne & Părău (Reference Guyenne and Părău2012, Reference Guyenne and Părău2014) for Toland’s model with reduced equations. The fundamental phenomenon of the nonlinear problem is the shedding of solitary waves. Here, we examine the problem by using the full Euler equations and Toland’s model. A single moving pressure which is defined as
$$\begin{eqnarray}P=A\text{e}^{-(x+220-ct)^{2}/16}\end{eqnarray}$$
is applied in (3.6) initially with amplitude
$A=0.2$
and translating speed
$c=1.3$
. The pressure is switched off at
$t=125$
. We found that the solution rapidly relaxes to a symmetric depression wave after the forcing is released (figure 14). The stability of this large-amplitude fully localised response is confirmed from its steady propagation at long times.
The stability analysis in § 5 shows that elevation waves with the structure of two well-separated troughs connected by a dimple are also longitudinally stable. We now show that these stable elevation waves can also be generated by moving loads. We choose the load as two fully localised pressures moving simultaneously on an ice sheet at the speed
$c$
. These pressures are defined as
$$\begin{eqnarray}\displaystyle & P_{1}=A\text{e}^{-(x+220-ct)^{2}/16}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & P_{2}=A\text{e}^{-(x+220-d-ct)^{2}/16}, & \displaystyle\end{eqnarray}$$
$A$
is a parameter measuring the magnitude of the pressures and
$d$
is the distance between the centres of the two pressures. Figure 15 shows snapshots of the generation of a stable elevation solitary wave by two pressures moving with the same speed. An elevation wave of two large troughs quickly forms after the forcing is turned off (figure 15
c,d). From longer computations (not shown), we believe that these solitary waves are stable and robust. This large response can propagate without changing its main structure, in spite of the interactions with the background radiation induced by the initial generation process. This numerical experiment can be used to mimic cars or aircraft running successively on a very large floating structure, and suggests that fully localised two-trough structure may be observable in a real situation.
Figure 15. Snapshots of the free surface due to two moving pressures (
$A=0.2$
and
$c=1.3$
) at
$t=25$
, 100, 140, 200, 325 (a–e). The forces are switched on at
$t=0$
and off at
$t=125$
. The centre of
$P_{1}$
is initially placed at
$x=-220$
and that of
$P_{2}$
is at
$x=-200$
.
7. Conclusions
Two-dimensional hydroelastic solitary waves propagating in deep water were investigated. The nonlinear model proposed by Toland (Reference Toland2008) was used to formulate the pressure exerted by the thin elastic sheet. There are depression and elevation solitary waves which are symmetric in the direction of wave propagation. These two branches, which were previously computed and extended by Guyenne & Părău (Reference Guyenne and Părău2012) and Wang et al. (Reference Wang, Vanden-Broeck and Milewski2013) respectively, were briefly reviewed. The depression and elevation branches bifurcate from the minimum of the phase speed and exist at subcritical speeds. The fact that different symmetric solitary waves can propagate at the same speed inspired us to search for asymmetric hydroelastic solitary waves. We superposed a depression wave and an elevation one as the initial guess to find an asymmetric solitary wave by using Newton’s method. The complete bifurcation diagram was then obtained by a numerical continuation method. These asymmetric waves have a two-packet structure and exist only in the subcritical regime. The asymmetric branch appears from a spontaneous symmetry-breaking bifurcation and ends with another symmetry-breaking bifurcation.
The stability properties of all of the hydroelastic solitary waves found to date, subject to longitudinal perturbations, were studied systematically by using the time-dependent conformal map technique together with a backward Euler’s method or a fourth-order Runge–Kutta method to integrate the full Euler equations in time. Our numerical experiments show that moderate-amplitude depression waves (lying between the stationary points (1) and (2) in figure 5) and certain elevation waves with two troughs connected by a dimple (lying in the segment between the turning point
$\unicode[STIX]{x2461}$
and the stationary point
$\unicode[STIX]{x2462}$
in figure 8) are stable. These two kinds of stable solitary waves can both be excited by applying one or two loads on the elastic sheet and moving the loads with a constant speed slightly below the phase speed minimum. Once the system gains enough localised energy, a stable hydroelastic solitary wave with a radiation field can be obtained after the pressure is switched off.
There are other possible fluid systems with complicated solitary-wave branches whose instability is worth investigating, most notably the interfacial gravity–capillary waves between two immiscible fluids. The results of Laget & Dias (Reference Laget and Dias1997) show that if the ratio of the density of the upper fluid to that of the lower fluid is larger than a critical value, the bifurcation diagram (figure 14 in their paper) is similar to figure 1. Therefore, a study of the stability of these interfacial solitary waves is of interest. It is also noted that the stability of multi-packet solitary waves was studied by Chardard, Dias & Bridges (Reference Chardard, Dias and Bridges2009, Reference Chardard, Dias and Bridges2011) using the Maslov theory. One might ask whether or not this theory can be applied to hydroelastic solitary waves.
We remark that we have considered in the present paper an idealised model which includes only the bending (i.e. beam-like response) of the material in order to identify the main properties of the waves. Other effects such as the stretching (i.e. rubber-like response), mass, thickness and inertia of the material can be considered in further studies. We finally comment that although overhanging waves are theoretically plausible in this idealised model, in practice the bending of the sheet will cause cracking at some critical curvature, and overturning of a finite-mass sheet must cause the sheet to fall off of the fluid.
Acknowledgement
This work was supported by EPSRC, under grant no. EP/J019569/1.
