]

# Nonlinear flexural-gravity waves due to a body submerged in the uniform
stream.

Y.A. Semenov  [  Institute of Hydromechanics of the National Academy of
Sciences of Ukraine, 8/4 Maria Kapnist Street, 03680 Kiev, Ukraine

###### Abstract

The two-dimensional nonlinear problem of steady flow past a body submerged
beneath an elastic sheet is considered. The mathematical model is based on the
velocity potential theory with fully nonlinear boundary conditions on the
fluid boundary and on the elastic sheet, which are coupled throughout the
numerical procedure. The integral hodograph method is employed to derive the
complex velocity potential of the flow which contains the velocity magnitude
on the interface in explicit form. The coupled problem has been reduced to a
system of nonlinear equations with respect to the unknown magnitude of the
velocity on the interface, which is solved using a collocation method. Case
studies are undertaken for both subcritical and supercritical flow regimes.
Results for interface shape, bending moment and pressure distribution are
presented for the wide ranges of Froude numbers and depths of submergence.
According to the dispersion equation, two waves on the interface may exist.
The first, longest wave is that caused by gravity, and the second, shorter
wave is that caused by the elastic sheet. The obtained solution exhibits
strongly nonlinear interaction of these waves above the submerged body. It is
found that near the critical Froude number, there is a range of submergences
in which the solution does not converge.

††preprint: AIP/123-QED

##  I INTRODUCTION

The problem of interaction between a fluid and an elastic boundary is a
classical problem of fluid mechanics which is of interest in offshore and
polar engineering, medicine and many industrial applications. In the last two
decades, this topic has received much attention due to the melting of ice in
the Arctic regions, opening new routes for ships and new regions for resource
exploration Squire et al. (1988, 1996); Korobkin et al. (2011). Most of the
studies are devoted to wave propagation along the ice sheet, its response on a
moving load, and effects of heterogeneous properties of the ice sheet such as
a floe, polynya, cracks, etc. Guyenne and Părău (2012, 2012, 2017); Li_2018 .

The works studying the interaction between the flow bounded by the ice sheet
and the body submerged in the fluid started relatively recently. Das and
Mandal Das_2006  considered oblique wave scattering by a circular cylinder
submerged beneath a uniform ice sheet of infinite extent and determined the
transition and reflection coefficients. To solve the problem, they employed
the multipole expansion method. Sturova Sturova_2015  applied the method of
matched eigenfunction expansions and studied the interaction of a submerged
cylinder and an inhomogeneous ice sheet including a floe or polynya. Tkacheva
Tkacheva_2015  considered oscillations of a cylindrical body submerged in a
fluid beneath an ice sheet and solved the problem through the Wiener-Hopf
technique. Savin and Savin Savin_2012  considered the ice perturbation by a
dipole submerged in water of infinite depth. They applied the complex variable
technique and the Fourier transform to solve the Laplace equation. Shishmarev
at al. Shishmarev_2019  studied the strains in an ice cover of a frozen
channel which are caused by a three-dimensional dipole moving under the ice at
a constant speed. Li et al. Li_2019  considered a circular cylinder submerged
below an ice sheet in water of finite depth. The solution method is based on
the derived Green function which satisfies the boundary conditions on the
ice/water interface. All the works mentioned above regarding submerged bodies
are based on linear potential flow theory, and the boundary value problem is
usually formulated in the frequency domain.

The nonlinear theory of hydroelasticity is currently under development. The
unknown shape of the ice/fluid interface and its higher-order derivatives,
which have to satisfy the dynamic boundary condition, are the main challenge
to derive analytical solutions or develop computational approaches. As the
dynamic boundary condition gets more complicates, e.g. include gravity,
surface tension and/or elasticity of the sheet covering the fluid, it
increases the level of the mathematical challenge which has to be addressed.
The simplest form of the dynamic boundary condition corresponds to free
streamline flows for which the velocity magnitude on the free streamline is
assumed to be constant. This class of flows is well developed and presented in
classical books by Milne-Thomson Milne-Thomson , Birkhoff and Zarantonello
Birkhoff , and Gurevich Gurevich .

For free-surface flows, gravity leads to an additional term in the dynamic
boundary condition which relates the velocity magnitude and the vertical
coordinate of the free surface. This kind of problem can be reduced to a
singular integrodifferential equation whose form depends on the solution
method and the choice of the governing functions. Various forms of the
integro-differential equation were derived by Forbes and Schwartz Forbes_1982
, King and Bloor King_Bloor , and Faltinsen and Semenov Faltinsen_Semenov .

For capillary free-surface flows, the dynamic boundary condition comprises the
curvature of the free surface which involves the first and second derivatives
of the free surface. Fewer analytical solutions for purely capillary waves are
presented in literature. Crapper Crapper  developed a closed-form solution for
a fluid of infinite depth, and Kinnersley Kinnersley  extended his method to a
fluid sheet. Crowdy Crowdy  developed a method based on complex function
theory and retrieved Kinnersley’s solution in much simpler form and obtained
new solutions for steady capillary waves on a fluid annulus. However, the
extension of the method to fluid/structure interaction problems with surface
tension seems to be nontrivial. Alternatively, several numerical methods have
been developed to solve the capillary and capillary-gravity flows. Schwartz
and Vanden-Broeck Schwartz_Vanden-Broeck  proposed a method based on a
boundary-integral formulation and finite difference approximation of the
derivatives and applied it to the purely capillary and capillary-gravity
waves. Vanden-Broeck and Miloh Vanden-Broeck_Miloh  proposed numerical methods
based on the Fourier-series expansion and studied steep gravity waves. Later,
the method was adopted by Blyth and Vanden-Broeck Blyth_Vanden-Broeck  and
Blyth and Părău Blyth_Parau  to compute nonlinear capillary waves on fluid
sheets of finite thickness. Yoon and Semenov Yoon_Semenov  considered a cavity
flow past a circular cylinder in the presence of surface tension and derived
the solution using the integral hodograph method. They derived a singular
integral equation with respect to the velocity magnitude, which is solved by
the method of successive approximations. The method can be applied to solve
problems with a more complicated form of the dynamic boundary condition which
comprises higher-order derivatives of the free surface. However, the higher-
order derivatives of the interface which appear in the dynamic boundary
condition results to a higher-order hypersingular integral equation. A special
numerical treatment is required to solve this type of integral equation.

The nonlinear theory of hydroelastic waves, for which the dynamic boundary
condition gets more complicated, has been studied intensively in recent
decades with emphasis on waves generated by a moving load. Most of the studies
are focused on the analysis and simulation of hydroelastic waves, which
account for the nonlinearity of the potential flow and elastic sheet
deformations. ParauPărău and Dias Parau_Dias  derived a forced nonlinear
Schrodinger equation for ice sheet deflection and studied the weakly nonlinear
effects. Fully nonlinear computations based on the boundary integral method
were presented by Guyenne and Părău Guyenne and Părău (2012). The nonlinear
response of an infinite ice sheet in the time domain has been studied by
Bonnefoy et al. Bonnefoy  using a higher-order spectral method. They found
that at a critical speed at which the linear response is infinite, the
nonlinear solution remains bounded. Despite the progress in the development of
numerical methods to solve nonlinear problems of hydroelastic waves, their
extension to a body submerged in fluid beneath an elastic sheet seems to be
not easy, since the flow potential has to satisfy an additional boundary
condition on the body surface.

In the present paper, we study a fully nonlinear problem of the hydroelastic
waves generated by a body submerged in the fluid beneath an elastic sheet. We
shall use the model of potential flow with the fully nonlinear kinematic and
dynamic boundary conditions on the submerged rigid body and the elastic sheet
which is modelled using the Cosserat theory of hyperelastic shells suggested
by Plotnikov and Toland Plotnikov . To solve the nonlinear problem, we adopt
the solution for a cylindrical body moving beneath a free surface Semenov_Wu ,
which is obtained using the integral hodograph method. An expression for the
flow potential which includes the velocity magnitude on the free surface in an
explicit form has been derived. This gives the possibility to adopt the
solution to in the present problem, because the velocity magnitude on the
interface between the fluid and elastic sheet appears in the dynamic boundary
condition explicitly. The coupling of the elastic sheet and fluid problems is
based on the condition of the same pressure on the interface, one from flow
dynamics and the second from elastic sheet equilibrium.

The derivation of the flow potential which contains in explicit form the
velocity magnitude on the interface and the numerical method to solve the
coupled fluid/elastic sheet interaction problem are presented in Section 2.
The extended numerical results are discussed in Section 3.

##  II THEORETICAL ANALYSIS

We consider a two-dimensional steady flow past a cylindrical body submerged
beneath an elastic sheet which is modelling the ice cover. The characteristic
length of the body is $L$, and the thickness of the sheet is $B_{i}$. A
Cartesian coordinate system $XY$ is defined with the origin at a point inside
the body and the $X$-axis along the velocity direction of the incoming flow
with a constant speed $U$. The $Y-$axis points vertically upwards. The fluid
is assumed to be inviscid and incompressible, and the flow is irrotational.
The elastic sheet is modelled by the Cosserat theory of hyperelastic shells
Plotnikov . The submerged rigid body is assumed to have an arbitrary shape
which can be defined by the slope of the body as a function of the arc length
coordinate $S$, or $\beta_{b}=\beta_{b}(S)$. The interface between the elastic
sheet and the liquid is defined by function $Y(X)$. The interactions of the
submerged body, flow and the elastic sheet may generate waves extending to
infinity in both upstream and downstream directions. On the other hand, the
flow is uniform at infinity, $Y\rightarrow\infty$ and $-\infty<X<\infty$, and
the velocity is constant there. In order to provide the same value of the flow
velocity at infinity in all directions, we introduce damping regions
$P_{1}P_{2}$ and $T_{1}T_{2}$ upstream and downstream, respectively, where a
term providing the wave damping is added in the dynamic boundary condition to
provide the same velocity magnitude at points $P_{2}$ and $T_{2}$, or
$V_{P2}=V_{T2}=U$. Outside the interval $P_{2}T_{2}$, the flow velocity on the
interface $V(X)\equiv U$ including infinity. Thus, the fluid surface has a
limit $Y(x)_{|x|\rightarrow\infty}=H$ which is defined as the submergence of
the cylinder measured from the origin of the coordinate system.

Figure 1:  ($a$) Physical plane and ($b$) the parameter $\zeta-$plane.

We will derive the complex potential of the flow, $W=W(Z)$, with $Z=X+iY$. For
the steady flow, the kinematic conditions on the body surface and the
interface mean that the stream function is constant, or
$\Im\\{W(Z)\\}=const.$, as they both are streamlines. The dynamic boundary
condition on the interface is obtained from the Bernoulli equation assuming
that the hydrodynamic pressure on the interface is the same as the pressure
conditioning the bending of the elastic sheet.

$\rho\frac{V^{2}}{2}+\rho gY+P_{ice}=\rho\frac{U^{2}}{2}+\rho gH+P_{a},$ (1)

where $U$ is the speed of the incoming flow, $\rho$ is the liquid density,
$V=|dW/dZ|$ is the magnitude of the complex velocity, $g$ is the acceleration
due to gravity, $H$ is the depth of submergence, and $P_{a}$ is the
atmospheric pressure. The pressure due to the bending of the elastic sheet is
Guyenne and Părău (2012); Plotnikov

$P_{ice}=D^{\prime}B^{3}_{i}\left(\frac{d^{2}\kappa}{dS^{2}}+\frac{1}{2}\kappa^{3}\right)+P_{a},$
(2)

where $D^{\prime}=E/(12(1-\nu))$ is the flexural rigidity of the elastic
sheet, $\kappa$ is the curvature of the interface, $B_{i}$ is the thickness of
the elastic sheet. Equation (2) corresponds to the assumptions: the elastic
sheet is inextensible and not prestressed Blyth_Parau .

Two different Froude numbers can be defined based on the characteristic length
$L$ or the depth of submergence $H$, respectively:

$F=\frac{U}{\sqrt{gL}},\qquad F_{h}=\frac{U}{\sqrt{gH}}.$ (3)

Using nondimensionalisation based on $U$, $L$ and $\rho$, we have $v=V/U$,
$x=X/L$, $y=Y/L$, $h=H/L$, $s=S/L$, $b_{i}=B_{i}/L$, and $W(Z)=ULw(z)$.
Replacing in equation (1) by (2), the dynamic boundary condition (1) takes the
form

$v^{2}=1-\frac{2(y-h)}{F^{2}}-2D\left(\frac{d^{2}\kappa}{dS^{2}}+\frac{1}{2}\kappa^{3}\right),$
(4)

where

$D=\frac{D^{\prime}}{\rho gLF^{2}},\qquad\kappa=\frac{d\delta}{ds}.$

The angle $\delta=\arcsin(dy/ds)=\pi+\beta$ is the angle between the $X-$axis
and the unit tangential vector $\mathbb{\tau}$ which is opposite to the
velocity direction $\beta$. The normal vector $\mathbb{n}$ is directed from
the liquid region outwards, while along the interface the spatial coordinate
increases in the direction of the vector $\mathbb{\tau}$ such that the liquid
region is on the left (see Fig. 1a).

Equation (4) contains the velocity magnitude along the interface, the wave
elevation $y$ and its derivatives which will be related in the following
throughout the derived expression for the flow potential.

###  II.1  Hodograph method.

Finding the function $w=w(z)$ directly is a complicated problem since the
boundary of the flow region is unknown. Instead, Joukovskii Joukovskii_1890
and Michell Michell  proposed to introduce an auxiliary parameter plane, or
$\zeta-$plane, which was typically chosen as the upper half-plane. Then, they
considered two functions, which were the complex potential $w$ and the
function $\omega=-\ln(dw/dz)$, both as functions of the parameter variable
$\zeta$. When $w(\zeta)$ and $\omega(\zeta)$ are derived, the velocity and the
flow region can be obtained in the parameter form as follows:

$\frac{dw}{dz}=\exp[-\omega(\zeta)],\qquad
z(\zeta)=z_{0}+\int_{0}^{\zeta}\frac{dw}{d\zeta^{\prime}}/\frac{dw}{dz}d\zeta^{\prime},$
(5)

where the function $z=z(zeta)$ is called as the mapping function.

The flow region beneath the interface and outside the body is a doubly
connected domain. A canonical region of a doubly connected domain is an
annulus. By making a cut connecting the external and the internal circles of
the annulus, the doubly connected region becomes simply connected. As shown in
Fig. 1a, $O_{-}D_{+}$ and $O_{+}D_{-}$ are the two sides of the cut which
could have an arbitrary shape but form a right angle with the flow boundary at
both the body surface (points $O_{-}$ and $O_{+}$) and the liquid/ice
interface (points $D_{-}$ and $D_{+}$).

The simply connected flow region $C_{-}D_{-}O_{+}BAO_{-}D_{+}C_{+}$ is then
transformed into the rectangular domain $O_{-}D_{+}CD_{-}BAO_{+}$ in the
parameter plane. An upper half-plane or unit circle is usually chosen as the
parameter plane. However, the flow region may have corner points at which the
mapping $z(\zeta)$ will be not conformal. Chaplygin (see chapter 1(5) in the
book Gurevich ) pointed out that there is flexibility in the choice of the
region of the parameter variable. It should be composed of straight lines and
arcs of circles in such a way that by means of image transformations of these
regions it is possible to cover the whole complex plane in a simple manner.

When solving boundary problems, the shape of the auxiliary parameter region is
usually chosen with the aim of obtaining the solution of a problem in the
simplest form with a minimal number of singular points at which the
transformation of the parameter region onto the complex potential region, $w$,
and the region of the function $dw/dz$ is not conformal. In the present case
of the doubly connected flow region, the additional corner points appear at
the intersections of two sides of the cut and the flow boundary. In order to
provide a conformal mapping at these corner points, we have chosen the
rectangle as the parameter domain, which also has right angles at points
$O_{-},O_{+},D_{-},D_{+}$, e.g. the same angles as in the physical plane.

When the parameter region is chosen as a half-plane or the first quadrant, the
polynomial functions are usually used to construct the mapping function
$z(\zeta)$. Here, for the rectangular domain, the polynomial functions will be
replaced by Jacobi’s theta functions Birkhoff , which are quasi-doubly-
periodic functions. Due the periodicity, they naturally satisfy the same
conditions on both sides of the cut. Jacobi’s functions have been used to
solve free surface problems involving doubly connected flow regions in the
books Milne-Thomson ; Birkhoff ; Gurevich ; Terentiev .

We may choose the coordinates of the rectangle vertices $O_{-}O_{+}D_{-}D_{+}$
as $(0,0)$, $(\pi,0)$, $(\pi,\pi\tau/4)$ and $(0,\pi\tau/4)$, respectively, as
shown in Fig. 1a. Here, $\tau$ is an imaginary number. The horizontal length
of the rectangle is then equal to $\pi$, and its vertical length is equal to
$\pi|\tau|/4$.

In the flow region, there are two stagnation points marked as $A$, where two
streamlines merge into one, and $B$, where a streamline splits into two
branches. Positions of these points in parameter plane $\zeta=a$, $\zeta=b$ as
well as the position of point $C$, $\zeta=c+\pi\tau/4$ which corresponds to
infinity in the physical plane. The parameters $a$, $b$ and $c$ should be
determined from additional conditions at the solution of the problem.

The interval $0\leq\xi\leq\pi$ on the real axis corresponds to the body
boundary. The interval $c<\xi\leq\pi$, $i\eta=\pi\tau/4$ corresponds to part
of the interface $D_{-}C_{-}$, and the interval $0\leq\xi<c$,
$i\eta=\pi\tau/4$ corresponds to the other part of the interface $D_{+}C_{+}$.
It should be noticed that points $C_{-}$ at $x\rightarrow-\infty$ and $C_{+}$
for $x\rightarrow+\infty$ in the physical plane have been transformed to the
same point $C$ in the parameter region $\zeta$.

###  II.2  Integral hodograph method:derivation of the governing functions
$dw/dz$ and $dw/d\zeta$.

At this stage we denote the angle of the velocity direction along the body as
$\beta_{b}(\xi)$ and the velocity magnitude on the free surface as $v(\xi)$.
With these notations, we have the following boundary-value problem for the
function of complex velocity, $dw/dz$:

$\left|\frac{dw}{dz}\right|=v(\xi),\qquad
0\leq\xi{\leq}\pi,\quad\eta=\pi\tau/4.$ (6)

$\chi(\xi)=\arg\left(\frac{dw}{dz}\right)=\left\\{{\begin{array}[]{l}-\beta_{b}(\xi),\qquad\quad
0\leq\xi<a,\quad\eta=0,\\\
-\beta_{b}(\xi)-\pi,\quad\,\,a<\xi<b,\quad\eta=0,\\\
-\beta_{b}(\xi)-2\pi,\quad b<\xi\leq\pi,\quad\eta=0.\end{array}}\right.$ (7)
$\frac{dw}{dz}(\xi=0,i\eta)=\frac{dw}{dz}(\xi=\pi,i\eta),\qquad 0\leq
i\eta\leq\pi\tau/4.$ (8)

In (7) the argument of complex velocity has the jumps equal to $-\pi$ at
stagnation points $A$ ($\zeta=a$) and $B$ ($\zeta=b$) due to the jump of the
velocity direction when passing through the stagnation point. The two vertical
sides of the rectangle in the parameter plane correspond to the two sides of
the cut in the physical plane. The velocities on both sides of the cut are the
same and therefore the condition of periodicity can be applied on the vertical
sides the rectangle. The solution of the boundary value problem (6)-(8) can be
obtained by applying the integral formulae derived in Sem_Wu2020 ,

$\displaystyle\frac{dw}{dz}$ $\displaystyle=$ $\displaystyle
v(\pi)\exp\left[-\frac{1}{\pi}\int_{0}^{\pi}\frac{d\chi}{d\xi}\ln\left(\frac{\vartheta_{1}(\zeta-\xi)}{\vartheta_{1}(\zeta-\xi-\pi\tau/2)}\right)d\xi\right.$
(9) $\displaystyle+$ $\displaystyle\frac{i}{\pi}\left.\int_{\pi}^{0}\frac{d\ln
v}{d\xi}\ln\left(\frac{\vartheta_{1}(\zeta-\xi-\pi\tau/4)}{\vartheta_{1}(\zeta-\xi+\pi\tau/4)}\right)d\xi+i\chi(\pi)\right].$

It can be easily verified that for $0<\xi<\pi$, $\eta=0$ the argument,
$\arg[(dw/dz)_{\zeta=\xi,\eta=0}]=\chi(\xi)$, while for $0<\xi<\pi$,
$i\eta=\pi\tau/4$, the modulus $|dw/dz|_{\zeta=\xi,i\eta=\pi\tau/4}=v(\xi)$,
i.e. the boundary conditions (6) and (7) are satisfied. The boundary condition
(8) is satisfied due to periodicity of the function $\vartheta_{1}(\zeta)$. By
substituting the boundary conditions (6) and (7) into (9) and evaluating the
first integral over the step change in the function $\chi(\xi)$ at points
$\zeta=a$ and $\zeta=b$, we obtain the expression for the complex velocity in
the rectangle $O_{-},O_{+},D_{=},D_{+}$ Sem_Wu2020 ,

$\displaystyle\frac{dw}{dz}$ $\displaystyle=$ $\displaystyle
v_{D}\frac{\vartheta_{1}(\zeta-a)\vartheta_{1}(\zeta-b)}{\vartheta_{4}(\zeta-a)\vartheta_{4}(\zeta-b)}\exp\left[\frac{1}{\pi}\int_{0}^{\pi}\frac{d\beta_{b}}{d\xi}\ln\frac{\vartheta_{1}(\zeta-\xi)}{\vartheta_{4}(\zeta-\xi)}d\xi\right.$
(10) $\displaystyle+$
$\displaystyle\frac{i}{\pi}\left.\int_{\pi}^{0}\frac{d\ln
v}{d\xi}\ln\frac{\vartheta_{1}(\zeta-\xi-\pi\tau/4)}{\vartheta_{4}(\zeta-\xi-\pi\tau/4)}d\xi-i\beta_{O}\right].$

where $\beta_{O}$ is the angle at point $O_{-}$ which is zero if point $O_{-}$
is the highest point of the body. The constant $v_{D}$ or the velocity
magnitude at point $D_{+}$, is determined by satisfying the velocity at
infinity, $\zeta=c+\pi\tau/4$, which is $1$, as it has been chosen as the
reference velocity, or

$\left|\frac{dw}{dz}\right|_{\zeta=c+\pi\tau/4}=1$ (11)

For steady flows, the stream function $\psi=\Im(w)$ takes constant values
along the body and the interface. According to Chaplygin’s special point
method Gurevich , to determine the function $w=w(\zeta)$, it is sufficient to
analyse all special points where the mapping is not conformal. These are the
stagnation points $A(\zeta=a)$ and $B(\zeta=b)$ and point
$C(\zeta=c+\pi\tau/4)$ corresponding to infinity in the $w-$plane. The order
of the function $w-w(\zeta)$ at these points can be determined by analysing
the behaviour of the argument of $w(\zeta)$ in the vicinity of these points.
Then, the derivation of the complex potential can be obtained in the form
Sem_Wu2020 .

$\displaystyle\frac{dw}{d\zeta}=K\frac{\vartheta_{1}(\zeta-a)\vartheta_{4}(\zeta-a)\vartheta_{1}(\zeta-b)\vartheta_{4}(\zeta-b)}{\vartheta_{1}^{2}(\zeta-c-\pi\tau/4)\vartheta_{1}^{2}(\zeta-c+\pi\tau/4)}.$
(12)

Dividing (12) by (10), we obtain the derivative of the mapping function as

$\displaystyle\frac{dz}{d\zeta}$ $\displaystyle=$
$\displaystyle\frac{K}{v_{D}}\frac{\theta_{4}^{2}(\zeta-a)\theta_{4}^{2}(\zeta-b)}{\theta_{1}^{2}(\zeta-c-\pi\tau/4)\theta_{1}^{2}(\zeta-c+\pi\tau/4)}$
$\displaystyle\times\exp\left[-\frac{1}{\pi}\int_{0}^{\pi}\frac{d\beta_{b}}{d\xi}\ln\frac{\theta_{1}(\zeta-\xi)}{\theta_{4}(\zeta-\xi)}d\xi\right.$
$\displaystyle-$ $\displaystyle\left.\frac{i}{\pi}\int_{\pi}^{0}\frac{d\ln
v}{d\xi}\ln\frac{\theta_{1}(\zeta-\xi-\pi\tau/4)}{\theta_{4}(\zeta-\xi-\pi\tau/4)}d\xi+i\beta_{O}\right].$

whose integration along the intervals $0\leq\xi<c$ and $c<\xi\leq\pi$ at
$\eta=\pi\tau/4$ in the $\zeta-$plane provides the parts $D_{+}C_{+}$ and
$D_{-}C_{-}$ of the free surface $C_{-}C_{+}$ in $\zeta-$plane, respectively.
The parameters $a,b,c,\tau$ and $K$, and the functions $\beta_{b}(\xi)$ and
$d(\ln v)/d\xi$ have to be determined from physical considerations and the
kinematic boundary condition on the body surface and the dynamic boundary
conditions on the free surface.

###  II.3  System of equations for parameters $a,b,c,\tau$ and $K$.

At infinity, point $C_{-}C_{+}(\zeta=c+\pi\tau/4)$, the velocity approaches
unit (since this velocity is chosen as the reference velocity), and its
direction is along the X-axis. Therefore, the argument of the complex velocity
(10) at point $\zeta_{C}=c+\pi\tau/4$ should be equal to zero

$\arg\left(\frac{dw}{dz}\right)_{\zeta=\zeta_{C}}=0.$ (14)

The scale factor $K$ is determined by the length $S_{b}$ which is the
perimeter of the body cross-section

$\int_{0}^{\pi}\frac{ds_{b}}{d\xi}d\xi=S_{b}.$ (15)

where

$\frac{ds_{b}}{d\xi}=\left|\frac{dz}{d\zeta}\right|_{\zeta=\xi}.$

The free surface on the left and right hand sides at infinity has the same
value of $y-$coordinate. This is also equivalent to that the stream function
$\psi=\Im(w)$ is continuous across the cut, or
$\Im(w_{D_{-}})-\Im(w_{D_{+}})=0$. By integrating $\Im(dw/d\zeta)$ along
$D_{-}D_{+}$ passing the point $\zeta_{C}$ along a semi-circle $C^{\prime}$ of
an infinitesimal radius $\varepsilon$, at which $dw/d\zeta$ in Eq.(12) has the
second order singularity, we have

$\displaystyle\Im\left(\int_{\pi}^{c+\varepsilon}\frac{dw}{d\zeta}d\zeta+\oint_{C^{\prime}}\frac{dw}{d\zeta}d\zeta+\int_{c-\varepsilon}^{0}\frac{dw}{d\zeta}d\zeta\right)$
$\displaystyle=$
$\displaystyle\Im\left(\oint_{C^{\prime}}\frac{dw}{d\zeta}d\zeta\right)=\Im\left(i\pi\begin{array}[]{c}{}{}{}\hfil\\\
\mbox{Res}\\\ ~{}^{\zeta=\zeta_{C}}\end{array}\frac{dw}{d\zeta}\right)$
$\displaystyle=$
$\displaystyle\Im\left\\{i\pi\frac{d}{d\zeta}\left[\frac{dw}{d\zeta}(\zeta-\zeta_{C})^{2}\right]_{\zeta=\zeta_{C}}\right\\}.$

Here the first and third terms on the left hand are zero because
$\Im(w)=const.$ on the free surface. From this equation it follows

$a+b=2c.$ (20)

The depth of submergence, $h$, and the flowrate, $Q$, between the body and the
free surface are related. Therefore, instead of a condition for the depth $h$,
we can use the following condition for the given flowrate $Q$, which is the
integral of the derivative of the complex potential along the side
$O_{-}D_{+}$

$\Im\left(\int_{0}^{\pi\tau/4}\frac{dw}{d\zeta}d\zeta\right)=Q.$ (21)

We may place a vortex with circulation $\Gamma$ at the centre of the cylinder,
which can be nondimensionalized as $\gamma=\Gamma/(2\pi UL)$. For a circular
cylinder, this does not affect the impermeable body surface boundary
condition, but does change the positions of the stagnation points and also
affects the free surface boundary. For a hydrofoil, $\gamma$ should be
determined through the Kutta condition at the trailing edge.

Integrating $dw/d\zeta$ along the body surface in the parameter plane, we have

$\Re\left(\int_{0}^{\pi}\frac{dw}{d\zeta}d\zeta\right)=2\pi\gamma.$ (22)

In the case $\gamma\neq 0$, the real part of the potential, $\phi=\Re(w)$,
have a jump on the sides $O_{-}D_{-}$ and $O_{+}D_{+}$ of the cut, while the
complex velocity, $dw/dz$ and the stream function $\psi=\Im(w)$ are still
continuous across the cut.

Equations (14) - (22) allow us to determine the unknown parameters
$a,b,c,\tau$ and $K$, which appear in the governing equations (10), (12) and
(II.2), once the functions $v(\xi)$ and $\beta_{b}(\xi)$ are specified.

###  II.4  Kinematic boundary conditions on the body surface.

By integrating the modulus of the derivative of the mapping function (II.2)
along the side $O_{-}O_{+}$ in the parameter plane, we can obtain the spatial
coordinate along the body as a function of the parameter variable

$s_{b}(\xi)=\int_{0}^{\xi}\frac{ds_{b}}{d\xi^{\prime}}d\xi^{\prime},$ (23)

where $ds_{b}/d\xi=|dz/d\zeta|_{\zeta=\xi,\eta=0}$. Since the function
$\beta_{b}(s_{b})$ is known, the function $\beta_{b}(\xi)$ can be determined
from the following integro-differential equation:

$\frac{d\beta_{b}}{d\xi}=\frac{d\beta_{b}}{ds_{b}}\frac{ds_{b}}{d\xi}.$ (24)

By substituting $dz/d\zeta$ from (II.2), this equation takes the form

$\displaystyle\frac{d\beta_{b}}{d\xi}$ $\displaystyle=$
$\displaystyle\kappa[s_{b}(\xi)]\frac{K}{v_{D}}\left|\frac{\theta_{4}^{2}(\xi-a)\theta_{4}^{2}(\xi-b)}{\theta_{1}^{2}(\xi-c-\pi\tau/4)\theta_{1}^{2}(\xi-c+\pi\tau/4)}\right|$
(25) $\displaystyle\times$
$\displaystyle\exp\left[-\frac{1}{\pi}\int_{0}^{\pi}\frac{d\beta_{b}}{d\xi^{\prime}}\ln\frac{\theta_{1}(\xi-\xi^{\prime})}{\theta_{4}(\xi-\xi^{\prime})}d\xi^{\prime}\right.$
$\displaystyle-$ $\displaystyle\left.\frac{i}{\pi}\int_{\pi}^{0}\frac{d\ln
v}{d\xi^{\prime}}\ln\frac{\theta_{1}(\xi-\xi^{\prime}-\pi\tau/4)}{\theta_{4}(\xi-\xi^{\prime}-\pi\tau/4)}d\xi^{\prime}\right],$

where $\kappa(s_{b})=d\beta_{b}/ds_{b}$ is the curvature of the body.

###  II.5  Nonlinear dynamic boundary condition.

The dynamic boundary condition (4) includes the interface elevation $y(s)$ and
its derivatives. Each branch of the interface, $C_{-}D_{-}(c<\xi<\pi)$ and
$C_{+}D_{+}(0<\xi<c)$, is evaluated by integration of the derivative of the
mapping function (II.2). The parameter form of the interface is as follows,

$s(\xi)_{\\{c<\xi\leq\pi,0\leq\xi<c\\}}=\int_{\\{\pi,0\\}}^{\xi}\frac{ds}{d\xi}d\xi,$
(26)
$y(\xi)_{\\{c<\xi\leq\pi,0\leq\xi<c\\}}=y_{D}+\Im\left(\int_{\\{\pi,0\\}}^{\xi}\left.\frac{dz}{d\zeta}\right|_{\zeta=\xi+\pi\tau/4}d\xi\right),$
(27)

where

$\displaystyle\frac{ds}{d\xi}$ $\displaystyle=$
$\displaystyle\left|\frac{dz}{d\zeta}\right|_{\zeta=\xi+\pi\tau/4}$ (28)
$\displaystyle=$
$\displaystyle\frac{K}{v(\xi)}\left|\frac{\vartheta_{4}^{2}(\xi-a+\pi\tau/4)\vartheta_{4}^{2}(\xi-b+\pi\tau/4)}{\vartheta_{1}^{2}(\xi-c)\vartheta_{4}^{2}(\xi-c)}\right|$

and $y_{D}$ is the vertical coordinate of points $D_{-}D_{+}$ and can be
obtained from

$y_{D}=\Im\left(i\int_{0}^{\pi|\tau|/4}\left.\frac{dz}{d\zeta}\right|_{\zeta=i\eta}d\eta\right).$
(29)

The curvature of the interface is

$\kappa[s(\xi)]=\frac{d\beta}{ds}=\frac{d\beta}{d\xi}/\frac{ds}{d\xi},$ (30)

where $d\beta/d\xi$ is determined by taking the argument of the complex
velocity from equation (10) at $\zeta=\xi+\pi\tau/4$,

$\displaystyle\beta(\xi)$ $\displaystyle=$
$\displaystyle\arg\left(\frac{dw}{dz}\right)=-\frac{1}{\pi}\int_{\pi}^{0}\frac{d\ln
v}{d\xi^{\prime}}\ln\left|\frac{\vartheta_{1}(\xi-\xi^{\prime})}{\vartheta_{4}(\xi-\xi^{\prime})}\right|-P_{1}(\xi),$
$\displaystyle P_{1}(\xi)$ $\displaystyle=$
$\displaystyle-\beta_{O}+\Im\left\\{\ln\frac{\vartheta_{1}(\xi-a+\pi\tau/4)\vartheta_{1}(\xi-b+\pi\tau/4)}{\vartheta_{4}(\xi-a+\pi\tau/4)\vartheta_{4}(\xi-b+\pi\tau/4)}\right\\}$
(32) $\displaystyle+$
$\displaystyle\frac{1}{\pi}\int_{0}^{\pi}\frac{d\beta_{b}}{d\xi^{\prime}}\Im\left\\{\ln\frac{\vartheta_{1}(\xi-\xi^{\prime}+\pi\tau/4)}{\vartheta_{4}(\xi-\xi^{\prime}+\pi\tau/4)}\right\\}d\xi^{\prime}$

and differentiating it respect to variable $\xi$,

$\frac{d\beta}{d\xi}=-\frac{1}{\pi}\int_{\pi}^{0}\frac{d\ln
v}{d\xi^{\prime}}\left(\frac{\vartheta_{1}^{\prime}(\xi-\xi^{\prime})}{\vartheta_{1}(\xi-\xi^{\prime})}-\frac{\vartheta_{4}^{\prime}(\xi-\xi^{\prime})}{\vartheta_{4}(\xi-\xi^{\prime})}\right)d\xi^{\prime}-P_{1}^{\prime}(\xi).$
(33)

Here, prime denotes the derivatives of the functions with respect to $\xi$.
The integrand of the above equation has a first-order singularity at point
$\xi^{\prime}=\xi$, since
$\vartheta_{1}(\xi-\xi^{\prime})\sim\xi-\xi^{\prime}$.

The derivatives of the curvature, $d\kappa/ds$ and $d^{2}\kappa/ds^{2}$, can
be obtained by differentiating (30). They will include higher-order
derivatives of the function $\beta(\xi)$ and a higher-order singularity in the
integrands, respectively. By substituting the derivatives of the curvature
into the dynamic boundary condition (4), we can derive a hypersingular
integral equation in terms of the function $d\ln v/d\xi$, the solution of
which requires special treatment. Instead of that, we will determine the
function $v(\xi)$ numerically using a collocation method.

###  II.6  Numerical method to determine the function $v(\xi)$.

If the tentative function $v(\xi)$ is given and the system of equations
(14)-(22) and the integrodifferential equation (25) are solved, then the
interface $z=z(\xi)$ depends only on the given function $v(\xi)$. We can chose
a fixed set of points $\hat{\xi}_{k}$, $k=1,\bar{K}$, distributed on the side
$D_{-}D_{+}$ of the parameter region corresponding to the interface. Then, the
nodes $s_{k}=s(\hat{\xi}_{k})$ and the coordinates $y_{k}=y(\hat{\xi}_{k})$ in
the physical plane are obtained using Eqs. (26) and (27). The curvature of the
interface and its derivatives are determined numerically applying spline
interpolation of the nodes $y_{k}$ in the intervals $(s_{k-1},s_{k})$. We
chose a fifth-order spline which provides continuous derivatives along the
interface up to the fourth derivative,

$\displaystyle y(s)=y_{k}+a_{1,k}(s-s_{k-1})+\ldots+a_{n,k}(s-s_{k-1})^{n},$
$\displaystyle\quad s_{k-1}<s<s_{k},\quad k=1,\ldots,\bar{K},\quad n=5.$ (34)

The curvature and its derivatives are obtained by differentiating Eq. (II.6):

$\beta=\arcsin
y^{\prime},\quad\kappa=\frac{y^{\prime\prime}}{\sqrt{1-y^{\prime
2}}},\quad\frac{d\kappa}{ds}=\frac{y^{\prime}y^{\prime\prime
2}-y^{\prime\prime\prime}(y^{\prime 2}-1)}{(1-y^{\prime 2})^{3/2}},\cdots.$

The system of nonlinear equations can be obtained by applying the dynamic
boundary condition (4) at the points $s_{k}=s(\hat{\xi_{k}})$, which is
written in the form

$\displaystyle G_{k}(\bar{V})=c_{pk}(\bar{V})-c_{pk}^{ice}(\bar{V})=0,\quad
k=1,\ldots,\bar{K},$ (35)

where $\bar{V}=(v_{1},\ldots,v_{\bar{K}})^{T}$ is the vector of unknown
velocities $v_{k}$ on the interface; the pressure coefficient on the interface
due to the flow is

$\displaystyle c_{pk}(\bar{V})=1-v_{k}^{2}-\frac{2[y(\bar{V})-h]}{F^{2}};$
(36)

and the pressure coefficient determining the bending of the elastic sheet is

$\displaystyle
c_{pk}^{ice}(\bar{V})=2D\left[\left(\frac{d^{2}\kappa}{ds^{2}}\right)_{k}+\frac{1}{2}\kappa_{k}^{3}\right].$
(37)

The system of equations (35) is solved using Newton’s method.The Jacobian of
the system is evaluated numerically using the central difference with $\Delta
v_{k}=10^{-5}$, $k=1,\bar{K}$. At each evaluation of the function
$G_{k}(\bar{V})$, the system of equations (14)-(22) and the
integrodifferential equation (25) are solved. From $5$ to $20$ iterations are
necessary to get convergence of the solution. All solutions, say
$\bar{V^{\ast}}$, reported here satisfied the condition

$\displaystyle\sum_{1}^{\bar{K}}|G_{k}(\bar{V^{\ast}})|<10^{-7}.$ (38)

which is regarded as giving a sufficiently accurate solution of the nonlinear
equations. However, the inaccuracy within the intervals
$(\hat{\xi}_{k-1},\hat{\xi}_{k})$ is somewhat lower. It will be discussed in
detail in section 3.1.

For supercritical flow regimes, the waving interface may extend to infinity.
However, due to the finite length of the calculation region and the condition
that the flow is uninform at infinity in all directions (upstream, downstream
and at infinite depth $y\rightarrow-\infty$), we need to introduce the damping
regions $P_{2}P_{1}$ upstream and $T_{1}T_{2}$ downstream. In these regions,
we add an artificial term in the boundary condition (4) which may be treated
as an external applied pressure,

$c_{p}=1-v^{2}-\frac{2(y-h)}{F^{2}}+C_{d}v\frac{dv}{ds},$ (39)

where the damping coefficient $C_{d}$ increases from $0$ at points $P_{1}$ and
$T_{1}$ to the values $C_{dL}$ and $C_{dR}$ at points $P_{2}$ and $T_{2}$,
respectively. The length of the damping regions are chosen to be about
$2\lambda_{0}$, where $\lambda_{0}$ is the wave length of the free surface
progressive waves according to the first approximation theory Lamb .

###  II.7  Dispersion equation.

Differentiating Eq. (4) with respect to the arc length coordinate along the
interface and taking into account that the slope of the interface
$\delta=\arcsin(dy/ds)$ is the angle between the unit tangential vector
$\mathbf{\tau}$ and the $x-$axis, we obtain

$F^{2}v^{2}\frac{d\ln
v}{ds}=-\sin\delta-F^{2}D\left[\frac{d^{4}\beta}{ds^{4}}+\frac{3}{2}\left(\frac{d\beta}{ds}^{2}\right)\frac{d^{2}\beta}{ds^{2}}\right].$
(40)

The above equation without an elastic sheet ($D=0$) and at small disturbances
of the free surface such that $\sin\delta\approx\delta$ can be written as

$v^{2}\frac{d\ln v}{ds}=-\frac{2\pi}{\lambda}\delta$ (41)

where the wave length $\lambda=\lambda_{0}=2\pi F^{2}$ is known from the
linear theory of gravity waves Lamb . In the presence of the sheet ($D\neq
0$), the waves of small amplitude approach a sinusoidal curve. Therefore,
their slope can be written as

$\delta(s)=\delta_{max}\cos\left(\frac{2\pi s}{\lambda}+\phi\right),$

where $\delta_{max}$ is the amplitude, and $\phi$ is the phase of the slope
relative to point $D_{+}$ at which $s=0$. Then, neglecting the square of the
curvature, i.e. the second term in brackets (40), we obtain

$v^{2}\frac{d\ln
v}{ds}=-\delta\left[\frac{1}{F^{2}}+D\left(\frac{2\pi}{\lambda}\right)^{4}\right].$
(42)

According to (41), the coefficient at the angle $\delta$ is equal to
$-2\pi/\lambda$, where $\lambda$ is the wave length of the interface.
Therefore, the following equation with respect to wave number $k=2\pi/\lambda$
is obtained

$k=\frac{1}{F^{2}}+Dk^{4}.$ (43)

This equation may have one, two or no roots, which depends on the constant $D$
depending on the thickness of the elastic sheet and the Froude number $F$. The
latter case corresponds to the subcritical flow regime for which the
perturbation on the interface decays in both directions. The case of one root
corresponds to the critical Froude number, $F_{cr}$, for which the waves of
the same length $\lambda_{cr}=2\pi/k_{cr}$ are extended to infinity in both
directions. Differentiating (43) with respect to $k$ and equating the result
to zero, after some manipulations we obtain

$k_{cr}=\sqrt[3]{4}\left(\frac{\rho
gL(1-\nu^{2})}{EB_{i}^{3}}\right)^{\frac{1}{4}},\quad
F_{cr}=\left(\frac{64}{81}\frac{EB_{i}}{\rho
gL(1-\nu^{2})}\right)^{\frac{1}{8}}.$ (44)

For $F>F_{cr}$, Eq. (43) has two roots, $k_{w}<k_{cr}$ due to gravity and
$k_{ice}>k_{cr}$ due to the elastic sheet. We note that the Eq. (42) is valid
along the whole interface, $-\infty<x<\infty$, so both waves associated with
wave numbers $k_{w}$ and $k_{ice}$ may appear upstream and downstream of the
submerged cylinder. We note that the depth of submergence of the body does not
appear in the dispersion equation, and therefore it does not influence the
wave number. However, we expect that the depth of submergence influences the
wave amplitude, similar to that observed for free-surface flows Sem_Wu2020 .

The roots of equation (43) for different Froude numbers and different
thicknesses of the elastic sheet are shown in figure 2.

Figure 2: Wave number vs. Froude number for different thicknesses of the
elastic: $b_{i}=0$ (solid line), $0.05$ (dashed line), $0.1$ (dotted line),
$0.2$ (dot-dashed line), $0.5$ (short dashed line).

Without an elastic sheet, the constant $D=0$, and Eq. (40) becomes $k=1/F^{2}$
that corresponds to the hyperbola in figure 2. For thickness $b_{i}>0$, the
critical Froude number appears as the minimal Froude to which corresponds only
one root and the vertical orientation of the slope. The larger relative
thickness of the elastic sheet results in the larger critical Froude number.
The dispersion equation predicts only possible waves, but the contribution of
each wave to the shape of the interface have to be determined from the
solution of the nonlinear problem.

##  III  Results

###  III.1  Numerical approach.

In the discrete form, the solution is sought on a fixed set of points
$\xi_{j}$, $j=1,N$, distributed along the side $O_{-}O_{+}$,
$0\leq\xi\leq\pi$, $\eta=0$, and a fixed set of points $\hat{\xi}_{i}$,
$i=1,M$ distributed along the side $D_{-}D_{+}$, $\eta=\pi/4$, of the
parameter region. The points $\xi_{j}$ are distributed so as to provide a
higher density of the points $s_{j}=s_{b}(\xi_{j})$ near stagnation points
$A(\zeta=a)$ and $B(\zeta=b)$. The distribution of the points $\hat{\xi}_{i}$
is chosen such to provide a higher density of the points
$s_{i}=s(\hat{\xi}_{i})$ closer to the body and the uniform distribution for
$|s_{i}|>\lambda$.

The number of nodes on the body and the interface are chosen in the ranges
$N=100\div 200$ and $M=1000\div 5000$, respectively, based on the requirement
to provide at least $80$ nodes within the shorter waves to get convergence and
reasonable accuracy of the solution. The nodes $\hat{\xi}_{k}$, $k=1,\bar{K}$,
used for interpolation of the interface, $y_{k}=y(\hat{\xi}_{k})$,
$s_{k}=s(\hat{\xi}_{k})$, are chosen on the set of points $\hat{\xi}_{i}$ such
that $i=4k$. Then, $\bar{K}=M/4=250\div 1000$ provides 20 nodes within the
shorter wave length at which the system of the nonlinear equations (35) is
solved.

The integrals appearing in Eq. (10) are evaluated based on the linear
interpolation of the functions $\beta_{b}(\xi)$ and $\ln v(\xi)$ on the
segments $(\xi_{j-1},\xi_{j})$ and $(\hat{\xi}_{i-1},\hat{\xi}_{i})$,
respectively. They are evaluated as follows:

$\displaystyle\frac{1}{\pi}\int_{0}^{\pi}\frac{d\beta_{b}}{d\xi}\ln\left(\frac{\vartheta_{1}(\zeta-\xi)}{\vartheta_{4}(\zeta-\xi)}\right)d\xi=\Delta\beta_{bj}[r_{j}(\zeta)+iq_{j}(\zeta)],$
$\displaystyle j=1,\ldots,N.\quad$ (45)
$\displaystyle\frac{i}{\pi}\int_{\pi}^{0}\frac{d\ln
v}{d\xi}\ln\left(\frac{\vartheta_{1}(\zeta-\xi-\pi\tau/4)}{\vartheta_{4}(\zeta-\xi-\pi\tau/4)}\right)d\xi=\Delta\ln
v_{i}[\hat{r}_{i}(\zeta)+i\hat{q}_{i}(\zeta)],$ $\displaystyle
i=1,\ldots,M.\qquad$ (46)

where $\Delta\beta_{bj}=\beta_{b}(\xi_{j})-\beta_{b}(\xi_{j-1})$, $\Delta\ln
v_{i}=\ln v(\hat{\xi}_{i})-\ln v(\hat{\xi}_{i-1})$,

$\displaystyle
r_{j}(\zeta)=\frac{1}{\pi\Delta\xi_{j}}\int_{\xi_{j-1}}^{\xi_{j}}\ln\left|\frac{\vartheta_{1}(\zeta-\xi)}{\vartheta_{4}(\zeta-\xi)}\right|d\xi,$
(47) $\displaystyle
q_{j}(\zeta)=\frac{1}{\pi\Delta\xi_{j}}\int_{\xi_{j-1}}^{\xi_{j}}\arg\left(\frac{\vartheta_{1}(\zeta-\xi)}{\vartheta_{4}(\zeta-\xi)}\right)d\xi,$
(48)
$\displaystyle\hat{r}_{i}(\zeta)=-\frac{1}{\pi\Delta\hat{\xi}_{i}}\int_{\hat{\xi}_{i-1}}^{\hat{\xi}_{i}}\arg\left(\frac{\vartheta_{1}(\zeta-\xi-\pi\tau/4)}{\vartheta_{4}(\zeta-\xi-\pi\tau/4)}\right)d\xi,$
(49)
$\displaystyle\hat{q}_{i}(\zeta)=\frac{1}{\pi\Delta\hat{\xi}_{i}}\int_{\hat{\xi}_{i-1}}^{\hat{\xi}_{i}}\ln\left|\frac{\vartheta_{1}(\zeta-\xi-\pi\tau/4)}{\vartheta_{4}(\zeta-\xi-\pi\tau/4)}\right|d\xi,$
(50)

The integrals (47) - (50) are evaluated using the $8-$point Legendre-Gauss
quadrature formula. The error of the solution within the intervals of
interpolation $(\hat{\xi}_{k-1},\hat{\xi}_{k})$ satisfies the relation

$\displaystyle\sum_{4k-4}^{4k}|G_{i}(\bar{V}^{\ast})|<10^{-3}$ (51)

which is regarded as giving a sufficiently accurate solution of the problem.

The method of successive approximations is adopted to solve the
integrodifferential equation (25), which in the discrete form becomes

$\frac{(\Delta\beta_{b})^{(k+1)}_{j}}{\Delta\xi_{j}}=\frac{\beta_{b}[s_{b}^{(k)}(\xi_{j})]-\beta_{b}[s_{b}^{(k)}(\xi_{j-1})]}{\Delta\xi_{j}},\qquad
j=1,\ldots,N,$ (52)

where the arc length along the body, $s_{b}^{(k)}(\xi)$, is evaluated using
(23) with $(\Delta\beta_{bj}/\Delta\xi_{j})^{(k)}$ known at the $(k)^{th}$
iteration. The iteration process converges very fast. After $5$ to $10$
iterations the error is below a prescribed tolerance of $10^{-6}$.

The derivative of the mapping function (II.2) has a second-order singularity
at point $\zeta=c+\pi\tau/4$. Therefore, points $\hat{\xi_{i}}$, $i=1,M$,
along the side $D_{-}D_{+}$ in the parameter region are distributed within the
two intervals $c+\varepsilon_{1}<\hat{\xi}_{i}\leq\pi$, $i=1,M_{1}$, and
$0\leq\hat{\xi}_{i}<c-\varepsilon_{2}$, $i=M_{1}+1,M$. These intervals
correspond to parts $D_{-}C_{-}$ and $D_{+}C_{+}$ of the interface
$C_{-}C_{+}$ in the physical plane. The values $\varepsilon_{1}$ and
$\varepsilon_{2}$ are chosen to provide the required length of the parts
$D_{-}C_{-}$ and $D_{+}C_{+}$.

###  III.2  Convergence study of the numerical method.

The formulation of the problem allows us to consider the free-surface flow
around the submerged circular cylinder if we chose zero thickness of the
elastic sheet. This case has been investigated in Sem_Wu2020  using the method
of successive approximations. The results based on the present collocation
method and that based on the successive approximations are shown in Fig. 3.

Figure 3: Verification of the numerical procedure comparing the shape of the
free surface using different methods: present collocation method (solid line);
successive approximation Sem_Wu2020  (dashed line); numerical solution Scullen
(symbols). Froude number $F=2.75$, depth of submergence $h=7.55$.

Without elastic sheet, the submerged body generates a progressive wave
downstream only. The free surface upstream for $x<-\lambda$ tends to be
parallel to the $x-$axis. This property was used in the model of Semenov & Wu
Sem_Wu2020  to determine a parameter which affects the velocity magnitude in
the damping region. Both the present numerical method and the computational
proceduresSem_Wu2020  predict the same shape of the free surface in the region
$s_{P1}<s<s_{T1}$. However, the present damping model provides the gradual
decay of the wave that is necessary to couple the pressure coefficients due to
the flow and the bending of the elastic sheet in the numerical procedure. It
is seen that the truncation of the computational region does not affect the
shape of the free surface in the interval $s_{P1}<s<s_{T1}$, as seen in Fig.
3.

Figure 4: Flexural-gravity waves generated by the submerged circular cylinder.
Thickness of elastic sheet $b_{i}=0.05$, Froude number $F=2$, critical Froude
number $F_{cr}=1.65$, depth of submergence $h=6.35$.

In order to investigate the effect of truncation in the presence of the sheet,
we consider an example of supercritical flow for the thickness $b_{i}=0.05$
and the depth of submergence $6.35$. From Eq. (43), we obtain the critical
Froude number, $F_{cr}=1.65$. For Froude number $F=2$, the flow is
supercritical, so there are two wave numbers $k_{w}=0.256$ and $k_{ice}=0.782$
determined from the dispersion equation (43). These wave numbers correspond to
the wave lengths $\lambda_{w}=24.5R$ and $\lambda_{ice}=8.03R$. The wave
length corresponding to the linear theory without an elastic sheet is
$\lambda_{0}=2\pi F^{2}R$. The interface is shown in Fig. 4 for two
computational regions $-8<x/\lambda_{0}<9$ (solid line) and
$-6<x/\lambda_{0}<7$ (dashed line). In the region $-5<x/\lambda_{0}<5$, where
the damping is absent ($C_{d}=0$) for both cases, the interfaces overlap.
Outside this region, the solid lines and the dashed lines start to diverge.
These results show that similar to free-surface flows, the truncation of the
computational region does not affect the part of the interface without
damping. The computations in Fig. 4 are carried out for the length of the
damping region, $x_{P1}-x_{P2}=\lambda_{0}$, and the damping region
downstream, $x_{T2}-x_{T1}=2\lambda_{0}$. The values $C_{dL}=2$ and
$C_{dR}=10$ were used, and the numbers of nodes are $N=100$ and $M=4000$.

###  III.3  Subcritical flows.

For Froude numbers $F<F_{cr}$, Eq. (43) has only complex roots, which
correspond to decaying perturbations of the interface caused by the submerged
cylinder. In Fig. 5, we show the interface profiles for the Froude number
$F=1.5$ and the relative thickness of the elastic sheet $b_{i}=0.05$ for
different depths of submergence. The interface shape is symmetric about the
$y-$axis. The shape is different from that observed for the free-surface flow
without an elastic sheet, for which the free surface is almost flat upstream
and exhibits a wave downstream. Thus, the elastic sheet supresses the waves
downstream and perturbs the flow upstream near the cylinder. As the thickness
of the elastic sheet tends to zero, the critical Froude number $F_{cr}$
decreases and become smaller than $F$. In this case, the flow becomes
supercritical, which drastically changes the interface shape. It will be
studied in the following subsection.

($a$)
($b$)

Figure 5: Free surface shape (a) and pressure coefficient (b) for Froude
number $F=1.5$, $F_{cr}=1.65$ and depth of submergence $h=5.37$ (solid line),
$5.60$ (dashed line), $6.35$ (dotted line) and $11.58$ (dash-dotted line).

It is expected that if cylinder is closer to the elastic sheet, or the depth
of submergence is smaller, then interaction between the cylinder and the
elastic sheet is stronger. This is observed in Fig. 5. The deflection of the
sheet above the cylinder exhibits a trough making the gap between the sheet
and the cylinder smaller. It is found that there is a minimal, or critical
depth of submergence, $h_{cr}$, below which the numerical solution cannot be
obtained. For $h$ slightly larger than $h_{cr}$, few iterations to solve the
system of nonlinear equations (II.6) are required to get the converged
solution, while for $h$ slightly smaller than $h_{cr}$, the elastic sheet
starts to oscillate, and the iterations do not converge.

Figure 6: Free surface shape (a) and pressure coefficient (b) for Froude
number $F=1.5$, $F_{cr}=1.65$ and depth of submergence $h=5.37$ (solid line),
$5.60$ (dashed line), $6.35$ (dotted line) and $11.58$ (dash-dotted line).

The interface profiles corresponding to for different Froude numbers are shown
in Fig. 6a. It is seen that the critical depth becomes larger as the Froude
number approaches the critical value . The corresponding distributions of the
pressure coefficient (the left axis) and the velocity magnitude (the right
axis) along the interface are shown in Fig. 6b. For the case of free-surface
flows, there is also a depth of submergence below which a steady free-surface
flow does not exist. In that case, the velocity magnitude at the crest of the
waves tends to zero, and the free surface shape forms a corner of . The
mechanism restricting the existence of the steady flow in the presence of the
elastic sheet is different. As shown in Fig. 6a, the velocity magnitude on the
interface is much larger than zero. The present results show that there is a
condition of consistency in interaction between the fluid, elastic sheet and
the submerged cylinder.

###  III.4  Supercritical flows.

We begin the computational analysis for relatively high Froude number, $F=3$,
or $F/F_{cr}=1.82$, and the ice thickness $b_{i}=0.05$, for which the
dispersion equation (43) has two real roots. The corresponding wave numbers
are $k_{ice}=1.13$ and $k_{w}=0.1113$. The second wave number almost coincides
with that obtained from the linear theory of gravity waves without an elastic
sheet, $k_{0}=1/F^{2}=0.1111$. The ratio of the wave lengths generated by the
submerged cylinder and the elastic sheet is $\lambda_{ice}/\lambda_{w}=10.11$.

The interfaces near the cylinder are shown in Fig. 7 for different depths of
submergence. The wave generated by the elastic sheet clearly seen upstream at
the smaller depths, and its amplitude decays as the depth of submergence
increases. For $h=6.13$ the wave upstream almost disappears. The wave
generated by the cylinder downstream is not completely seen since its length
exceeds the length of the interface shown in Fig. 7. The larger length of the
interface is shown in Fig. 8. For submergence $h=6.13$, the interface
coincides with that obtained without an elastic sheet which is shown in Fig. 8
by the dashed line.

Figure 7: Interface shape at different depths of submergence for Froude number
$F=3$ and thickness of elastic sheet $b_{i}=0.05$.

As the submergence decreases, the amplitude of the wave downstream increases
and reaches its maximal value at $h=2.95$, and then it starts to decrease.
This feature has been studied in Sem_Wu2020 . As the cylinder approaches the
free surface, it affects the free surface at smaller distances from the
cylinder, and the flow tends to be symmetric about the $y-$axis, similar to
that for $F\rightarrow\infty$.

The elastic sheet weakly influences the interface downstream but generates
wave upstream. As the cylinder approaches the free surface, the bending of the
elastic sheet near the cylinder increases, as can be seen from Figs. 7 and 8.
This causes the increase of the amplitude of the wave upstream.

Figure 8: Effect of submergence on the interface shape for the circular
cylinder at Froude number $F=3$. The thickness of elastic sheet $b_{i}=0.05$
(solid lines) and $b_{i}=0$ (dashed lines).

The bending moment and the pressure coefficient are shown in Figs. 9 and 10
for depths of submergence and , respectively. Although the wave due to the
elastic sheet is invisible (Fig. 9$c$), the pressure coefficient and the
bending moment oscillate at both directions upstream (left axis in Fig. 9$b$)
and downstream of the cylinder (right axis in Fig. 9$b$). The frequency of
oscillations at the upstream is caused by the elastic sheet, while the
frequency of oscillations at the downstream is due to gravity. This
qualitatively agrees with results based on the linear theory Savin_2012 ;
Li_2019 .

Figure 9: Bending moment ($a$), pressure coefficient ($b$) and interface shape
($c$) at Froude number $F=3.0$. The thickness of elastic sheet $b_{i}=0.05$
and depth of submergence $h=6.13$.
Figure 10: The same as in Fig. 9 at the depth of submergence $h=1.94$.

The results for Froude number $F=2.5$, or $F/F_{cr}=1.52$, and ice thickness
$b_{i}=0.05$ are shown in Figs. 11-14. The wave numbers are as follows:
$k_{ice}=0.97$, $k_{w}=0.1606$ and $k_{0}=0.1600$. The ratio of the wave
lengths $\lambda_{w}/\lambda_{ice}=6.04$. The interface shapes for different
depths of submergence are shown in Fig. 11 and 12. They are similar to those
in Figs. 7 and 8 for $F=3.0$. However, the pressure coefficient and the
bending moment shown in Fig. 13$a$,$b$ exhibit behaviour corresponding to
superposition of the gravity wave (longer wave) and the elastic sheet wave
(shorter wave). The amplitude of the bending moment corresponding to the
gravity wave is larger than that corresponding to the elastic sheet wave. The
latter becomes largest at the trough of the gravity wave, and it almost
disappears at the crest, as seen in Figs. 13$b$ and 14$b$. Such behaviour of
the bending moment and the pressure coefficient demonstrates the nonlinear
interaction of the elastic sheet and the flow, which is still invisible for
the interface profile in Figs. 13$c$ and 14$c$.

Figure 11: Interface shape at different depths of submergence for Froude
number $F=2.5$ and thickness of elastic sheet $b_{i}=0.05$.
Figure 12: Effect of submergence on the interface shape for the circular
cylinder at Froude number $F=2.5$.
Figure 13: Bending moment ($a$), pressure coefficient ($b$) and interface
shape ($c$) at Froude number $F=2.5$ and the submergence $h=4.4$.
Figure 14: The same as in Fig. 13 at the depth of submergence $h=1.94$.

For Froude number $F=2$, or $F/F_{cr}=1.21$, the wave numbers are as follows:
$k_{ice}=0.782$, $k_{w}=0.256$ and $k_{0}=0.250$. The ratio of the wave
lengths $\lambda_{w}/\lambda_{ice}=3.05$. The interface shapes for depths of
submergence in the range from $4.6$ to $11.6$ are shown in Fig. 15 and 16.

At the upstream direction, the wave caused by the elastic sheet becomes
visible even for the relatively large depth of submergence, $h=11.6$. Its
shape is like a sinusoid with constant amplitude. The amplitude increases as
the depth of submergence decreases. At depth of submergence $h=4.61$, the
elastic sheet interacts with the flow in such a way that the wave due to
gravity extends in the upstream direction, and the interface exhibits
superposition of the both waves.

At the downstream direction, the interface shape differs from a wave of
constant amplitude. The shape near the body corresponds to the superposition
of the gravity and elastic sheet waves, and then gradually approaches to the
pure gravity wave. This also occurs for larger submergence but is less
visible. The contribution of the elastic sheet wave decays downstream, and the
interface approaches the pure gravity wave far downstream.

Figure 15: Interface shape at different depths of submergence for Froude
number $F=2.0$.
Figure 16: Effect of submergence on the interface shape for the circular
cylinder at Froude number $F=2.0$.
Figure 17: Bending moment ($a$), pressure coefficient ($b$) and interface
shape ($c$) at Froude number $F=2.0$ and the submergence $h=11.6$.
Figure 18: The same as in Fig. 17 at the depth of submergence $h=4.61$.

The pressure coefficient and the bending moment along the interface are shown
in Fig. 17$a,b$ for $h=11.65$ and in Fig. 18$a,b$ for $h=4.6$. They
demonstrate interaction between the gravity and the elastic sheet waves. The
wave due to the elastic sheet is dominating near the cylinder. It keeps
constant amplitude in the upstream direction and gradually decays in the
downstream direction. For $x/\lambda_{0}>8$, the amplitude of the bending
moment in Fig. 18$a$ caused by the elastic sheet becomes smaller than that
caused by the gravity wave, and the oscillations become qualitatively similar
to that in Figs. 13$a$ and 14$a$. It is seen from Fig. 18$a$ that the distance
at which the bending moment caused by the elastic sheet decays is much larger
than that for $F=2.5$.

Figure 19: Interface shape at different depths of submergence for Froude
number $F=1.7$.
Figure 20: Effect of submergence on the interface shape for the circular
cylinder at Froude number $F=1.7$.
Figure 21: Bending moment ($a$), pressure coefficient ($b$) and interface
shape ($c$) at Froude number $F=1.7$ and the submergence $h=11.6$.
Figure 22: The same as in Fig. 21 at the depth of submergence $h=6.8$.

The results for Froude number $F=1.7$, which is quite close to the critical
Froude number, $F/F_{cr}=1.03$, are shown in Figs. 19-23. The wave numbers
$k_{ice}$ and $k_{w}$ approach each other, and the ratio of the wave lengths
becomes smaller, $\lambda_{w}/\lambda_{ice}$. The interfaces are shown in
Figs. 19 and 20 for the depths of submergence from $11.6$ to the $6.8$. By
comparing the interfaces at the depth $h=11.6$ for the Froude numbers $2.5$
and $2$ in Figs. 12 and 16, respectively, we can find that the amplitude of
the interface wave increases at the region upstream while the amplitude at the
region downstream of the cylinder decreases. The latter agrees with the free-
surface gravity flow past the submerged circular cylinder Sem_Wu2020 . The
former indicates the larger effect of the elastic sheet at the smaller Froude
number.

The shapes of the interface for different depths of submergence in Fig. 20 are
similar each other, but the amplitudes are different. In the upstream
direction, $x/\lambda_{0}<0$, the shapes are periodic with period of about
$2\lambda_{0}$, that corresponds to the superposition of two sinusoidal waves
with the gravity and elastic sheet waves.

In the downstream direction, $x/\lambda_{0}>0$, the shape of the interface is
not exactly periodic because the amplitude of the elastic sheet wave decays
slowly downstream. By comparing the shape of the interfaces, behaviour of the
bending moment and the pressure coefficient for different Froude numbers, we
can see that the wave due to the elastic sheet decays slower as the Froude
number approaches the critical value. The similar behaviour of the bending
moment and the pressure coefficient can be seen in Figs. 21 and 23.

Figure 23: Onset of the solution convergence.

The interaction between the fluid and elastic sheet and between the fluid and
the submerged body is not always consistent for the steady flow. This
situation has some analogy with the free-surface flow past the submerged
cylinder, for which there is some combination of the Froude number and depth
of submergence at which the steady solution does not exists. As the depth of
submergence increases, the region of non-existence of the solution decreases
and then disappears at large enough submergence. There is a maximal possible
deflection of the free surface at which the dynamic boundary condition can be
satisfied. In the presence of the elastic sheet, the situation becomes more
complicated because the dynamic boundary condition includes not only the
deflection of the interface but also its derivatives. The term due to gravity
in Eq. (31) increases as the Froude number decreases. It may cause such
combination of the deflection, its derivatives and the velocity magnitude on
the interface that the pressure distribution for the fluid (31) and for the
elastic sheet (32) cannot be the same, or the dynamic boundary condition (4)
cannot be satisfied. For a supercritical flow, the fluid forces dominate the
elastic forces, while for subcritical flows vice versa. This changes the flow
configuration near the critical Froude number. The interface shown in Fig. 19
for F = 1.7 and corresponds to the supercritical flow closest to the onset
where the converged solution can be obtained. It forms a hill over the
cylinder, while for the subcritical flow in Fig. 6 for , the interface forms a
trough. The different limits of the flow configurations for sub- and
supercritical regimes indicate an inconsistency of the fluid and elastic
forces in some range near the critical Froude number. As the submergence
increases, the deflection of the interface decreases as well as the range of
Froude numbers at which the steady flow and the elastic sheet are not
consistent. That can be seen in Fig. 23 where the onset of the region of
existence of the steady solution is shown in the plane of the parameters
Froude number vs. depth of submergence.

##  IV  Conclusions

A fully nonlinear numerical solution for the problem of steady gravity flow
past a body submerged beneath an elastic sheet is presented in the form of the
nonlinear analytical solution for the fluid part of the problem and the
nonlinear Cosserat plate model applied to the elastic sheet, which are coupled
throughout the numerical procedure. The solution of the fluid part of the
problem is based on the integral hodograph method employed to construct the
complex potential of the flow and Jacobi’s elliptic theta functions to deal
with the doubly connected fluid domain. The curvature and higher-order
derivatives of the fluid boundary involved in the nonlinear Cosserat plate
model have been evaluated using spline interpolation. The coupled problem has
been reduced to a system of nonlinear equations with respect to the unknown
magnitude of the velocity on the interface, which are solved using a
collocation method.

Steady solutions of the full nonlinear problem were computed for sub- and
supercritical regimes. For subcritical regimes, the dispersion equations have
no real roots. The elastic sheet exhibits a most considerable deflection above
the cylinder, which rapidly decays away. The deflection forms a curve
symmetric about the Y-axis with a trough above the cylinder. The trough
becomes larger as the depth of submergence decreases. These results
qualitatively agree with linear solutions Savin_2012 ; Li_2019 . At the same
time, the present nonlinear solution revealed a critical depth below which the
deflection of the interface cannot provide balance between the bending and
hydrodynamic forces in the steady flow. The critical submergence increases as
the Froude number approaches the critical Froude number.

For supercritical regimes, the dispersion equations have two positive real
roots which correspond to two wave numbers. The smaller wave number is closer
to that corresponding to the gravity wave behind the cylinder without an
elastic sheet, and the larger wave number appears in the presence of the
elastic sheet. The dispersion equation does not restrict the flow regions in
which the waves may occur, i.e. upstream or downstream of the cylinder. From
the linear theories Savin_2012 ; Li_2019 , it was found that the gravity wave
occurs downstream of the cylinder, while the wave due to the elastic sheet
occurs upstream. The present nonlinear solution revealed that the waves may
occur in both directions, but their amplitudes in each direction significantly
depend on the perturbation of the interface and the ratio $F/F_{cr}$.

The calculations are presented for the thickness of the elastic sheet
$b_{i}=0.05$ that corresponds to the critical Froude number $F_{cr}=1.65$. For
larger thickness of the sheet, the critical Froude number increases. We expect
that the flow configurations for larger thicknesses of the ice sheet will be
similar to those for $b_{i}=0.05$ at the same ratio $F/F_{cr}$.

At high Froude number $F/F_{cr}>1.8$ and depth of submergence $h>6$, the
interface shape is almost the same as it is without the elastic sheet.
However, the effect of the elastic sheet can be seen in the behaviour of the
bending moment and the pressure coefficient at the upstream. They oscillate
with the wave number corresponding to the elastic sheet. The amplitude of
oscillations is larger than that corresponding to the gravity wave at the
downstream. At smaller submergence, the perturbation of the interface
increases, and it becomes visible at the upstream. At very small depth of
submergence, the elastic sheet starts to affect the whole interface: at the
downstream, the amplitude of the gravity wave slightly decreases; and at the
upstream, the gravity wave is excited in addition to the elastic sheet wave.
The interface represents superposition of both these waves.

As the Froude number decreases and approaches the critical value $F_{cr}$, the
wave caused by the elastic sheet can be observed in both directions. Its
amplitude is constant in the upstream direction. In the downstream direction,
the contribution of the elastic wave to the resulting shape decays. The closer
the Froude number $F$ to the critical value $F_{cr}$, the slower decay is
observed.

Similar to subcritical regimes, there is a critical submergence below which
the steady supercritical solution cannot be obtained. The closer the Froude
number to the critical value $F_{cr}$, the larger the critical depth of
submergence. This may be caused due to the interaction of the two waves with
closer wave length that may cause a resonance phenomenon.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Squire et al. (1988) V. A. Squire, W. H. Robinson, P. J. Langhorne, and T. G. Haskell, "Vehicles and aircraft on floating ice," Nature 333, 159-161 (1988).
  * Squire et al. (1996) V. A. Squire, R. J. Hosking, A. D. Kerr, and P. J. Langhorne, "Moving Loads on Ice Plates", Kluwer, 1996.
  * Korobkin et al. (2011) A. Korobkin, E. I. Părău, and J.-M. Vanden-Broeck, "The mathematical challenges and modelling of hydroelasticity," Phil. Trans. R. Soc. Lond. A 369, 2803-2812 (2011).
  * Guyenne and Părău (2012) P. Guyenne and E. I. Părău, "Computations of fully nonlinear hydroelastic solitary waves on deep water," J. Fluid Mech. 713, 307-329 (2012).
  * Guyenne and Părău (2017) P. Guyenne and E. I. Părău, "Numerical study of solitary wave attenuation in a fragmented ice sheet," Phys. Rev. Fluids, 2, 034002 (2017).
  * (6) Z. F. Li, G. X. Wu and C. Y. Ji, "Wave radiation and diffraction by a circular cylinder submerged below an ice sheet with a crack," J. Fluid Mech. 845, 682-712 (2018).
  * (7) D. Das and B. N. Mandal, "Oblique wave scattering by a circular cylinder submerged beneath an ice-cover," Int. J. Eng. Sci. 44, 166-179 (2006).
  * (8) I. V. Sturova, "Radiation of waves by a cylinder submerged in water with ice floe or polynya," J. Fluid Mech. 784, 373-395 (2015).
  * (9) L. A. Tkacheva, "Oscillations of a cylindrical body submerged in a fluid with ice cover," J. Appl. Mech. Tech. Phys. 56, 1084-1095 (2015).
  * (10) A. A. Savin and A. S. Savin, "Ice cover perturbation by a dipole in motion within a liquid," Fluid Dyn. 47, 139-146 (2012).
  * (11) K. Shishmarev, T. Khabakhpasheva, A. Korobkin, "Ice response to an underwater body moving in a frozen channel," Applied Ocean Research 91, 101877 (2019).
  * (12) Z. F. Li, G. X. Wu, and Y. Y. Shi, "Interaction of uniform current with a circular cylinder submerged below an ice sheet," Appl. Ocean Res. 86, 310-319 (2019).
  * (13) L. M. Milne-Thomson, "Theoretical Hydrodynamics," 5th Edition (Dover Publications, 1968).
  * (14) G. Birkhoff and E. H. Zarantonello, "Jets, Wakes and Cavities", (Academic Press, 1957).
  * (15) M. I. Gurevich, Theory of Jets in Ideal Fluids (Academic Press, 1965).
  * (16) L. K. Forbes and L. W. Schwartz, "Free-surface flow over a semicircular obstruction," J. Fluid Mech. 114, 299-314 (1982).
  * (17) A. C. King and M. I. G. Bloor, "A semi-inverse method for free-surface flow over a submerged body," Q. J. Mech. Appl. Maths 42, 183-202 (1989).
  * (18) O. M. Faltinsen and Y.A. Semenov, "The effect of gravity and cavitation on a hydrofoil near the free surface," J. Fluid Mech. 597, 371-394 (2008).
  * (19) G. D. Crapper, "An exact solution for progressive capillary waves of arbitrary amplitude," J. Fluid Mech. 2, 532-540 (1957).
  * (20) W. Kinnersley, "Exact large amplitude capillary waves on sheets of fluid," J. Fluid Mech., 77, 229-241 (1976).
  * (21) D. G. Crowdy, "Exact solutions for steady capillary waves on a fluid annulus," J. Nonlinear Sci. 9, 615-640 (1999).
  * (22) L. W. Schwartz and J.-M. Vanden-Broeck, "Numerical solution of the exact equations for capillary-gravity waves," J. Flud Mech. 95, 119-139 (1979).
  * (23) J.-M. Vanden-Broeck and T. Miloh, "Computations of steep gravity waves by a refinement of Davies-Tulin’s approximation," SIAM J. Appl. Math. 55 (4), 892-903 (1995).
  * (24) M. G. Blyth and J.-M. Vanden-Broeck, "New solutions for capillary waves on fluid sheets," J. Fluid Mech. 507, 255-264 (2004).
  * (25) M. G. Blyth, E. I. Părău, and J.-M. Vanden-Broeck, "Hydroelastic waves on fluid sheets," J. Fluid Mech. 689, 541-551 (2011).
  * (26) B-S. Yoon and Y. A. Semenov, "Capillary cavity flow past a circular cylinder," Eur. J. Mech. B Fluids 28, 670-676 (2009).
  * (27) E. I. Părău and F. Dias, "Nonlinear effects in the response of a floating ice plate to a moving load, "J. Fluid Mech. 460, 281-305 (2002).
  * (28) F. Bonnefoy, M. H. Meylan, and P. Ferrant, Nonlinear higher-order spectral solution for a two-dimensional moving load on ice, J. Fluid Mech. 621, 215-242 (2009).
  * (29) P. I. Plotnikov and J. F. Toland, "Modelling nonlinear hydroelastic waves," Phil. Trans. R. Soc. Lond. A 369, 2942-2956 (2011).
  * (30) Y. A. Semenov and G. X. Wu, "Free-surface gravity flow due to a submerged body in uniform current," J. Fluid Mech. 883, A60 (2020).
  * (31) N. E. Joukovskii, "Modification of Kirchhoff’s method for determination of a fluid motion in two directions at a fixed velocity given on the unknown streamline", Math. Sbornik. 15 (1), 121-278 (1890).
  * (32) J. H. Michell, "On the theory of free stream lines", Phil. Trans. R. Soc. Lond. A 181, 389-431 (1890).
  * (33) J. S. Terentiev, A. G. Kirschner, and I. N. Uhlman, "The Hydrodynamics of Cavitating Flow." (Backbone Publishing Company, 2011).
  * (34) Y. A. Semenov and G. X. Wu, "Free-surface gravity flow due to a submerged body in uniform current", J. Fluid Mech. 883, A60 (2020).
  * (35) H. Lamb, _Hydrodynamics._ Cambridge University Press, (1932).
  * (36) D. Scullen and E. O. Tuck, "Nonlinear free-surface flow computations for submerged cylinders", J. Ship Res. 39, 185-193 (1995).
  * Smith et al. (2018) F. Smith, A. Korobkin, E. Părău, D. Feltham, and V. Squire, "Modelling of sea-ice phenomena," Phil. Trans. R. Soc. Lond. A 376 , 20180157 (2018).