1. Introduction
In recent years, new commercial shipping routes through the Arctic Ocean have gradually become visible (Smith & Stephenson Reference Smith and Stephenson2013). Icebreakers are commonly employed to ensure safe and environmentally friendly navigation for offshore operations. However, conventional icebreakers are often unsuitable for narrow channels, where limited manoeuvring space presents operational challenges, and are inefficient in thin ice conditions due to high fuel consumption and environmental impact. A promising alternative involves utilising flexural gravity waves generated by submerged vehicles (Kozin & Onishchuk Reference Kozin and Onishchuk1994) or moving hovercraft on the ice sheet (Eyre Reference Eyre1977). To apply these methods effectively, it is crucial to understand the physical mechanisms of the hydrodynamic interactions between the moving body, flexural gravity waves and ice sheet deformation.
For loads moving on an ice sheet, Takizawa (Reference Takizawa1985) conducted an experiment on an ice-covered lake, and illustrated that the distribution of the generated flexural gravity waves in the ice sheet was related to the speed of the load. To explain this observed phenomenon, the ice sheet may be modelled as a Kirchhoff–Love plate, and the linearised velocity potential theory for the fluid may be adopted. Based on this model, Davys, Hosking & Sneyd (Reference Davys, Hosking and Sneyd1985) employed the Fourier transform method and studied the three-dimensional (3-D) problem of hydroelastic waves generated by a single source moving steadily on an ice sheet. Later, Schulkes & Sneyd (Reference Schulkes and Sneyd1988) studied the 2-D transient problem of waves generated by moving loads. Their results showed that, for this type of problem, there are two critical Froude numbers
$F=F_c$
and the depth-based
$F=1$
. When
$F\lt F_c$
, there is no travelling wave to infinity and only evanescent ones; when
$F_c\lt F\lt 1$
, two travelling waves emerge, travelling upstream and downstream, respectively; when
$F\gt 1$
, only the upstream wave remains. Furthermore, Milinazzo, Shinbrot & Evans (Reference Milinazzo, Shinbrot and Evans1995) and Nugroho et al. (Reference Nugroho, Wang, Hosking and Milinazzo1999) considered the problem of flexural waves generated by prescribed point and distributed loads at the critical Froude numbers. They found that the velocity potential and ice sheet deflection for both the 2-D and 3-D problems are unbounded at
$F=F_c$
. In addition to a steadily moving load, Miles & Sneyd (Reference Miles and Sneyd2003) investigated the response of a floating ice sheet to a line load accelerating from rest. Their results indicated that the hydroelastic wave profile underwent a stable transition when the speed of the load passed through the critical speed
$F=F_c$
. More recently, Hosking & Milinazzo (Reference Hosking and Milinazzo2022) further extended the work to consider line loads of arbitrarily varying speed.
For bodies moving below an ice sheet, Savin & Savin (Reference Savin and Savin2012) simplified a 2-D circular cylinder as a dipole, and obtained an analytical solution for a dipole moving uniformly beneath a homogeneous ice sheet. Sturova (Reference Sturova2013) further solved the 3-D problem of a sphere moving at a forward speed based on the multipole expansion procedure (Wu Reference Wu1995). Li, Wu & Shi (Reference Li, Wu and Shi2019) considered an equivalent problem of a uniform current interaction with a circular cylinder submerged below an ice sheet. The Green function or the velocity potential due to a single source was first derived, and then multipole expansion was adopted to construct the velocity potential. They noticed that the resistance and lift on the body experienced very rapid change when
$F$
was near
$F_c$
or
$1$
. These were quite different from the free surface problem, where there is only one critical speed at
$F=1$
(Lighthill Reference Lighthill1978). Later, Yang, Wu & Ren (Reference Yang, Wu and Ren2021) further extended the work to the 3-D ice-covered channel problem and found that the Green function is singular at an infinite number of critical speeds. A major problem when the Green function becomes singular is that the commonly used boundary element method (BEM) is no longer solvable. The behaviour of the flow corresponding to the real body then becomes unknown because the solution cannot be found. We may notice that the present case is very much different from that of a load moving on the ice sheet (Milinazzo et al. Reference Milinazzo, Shinbrot and Evans1995; Nugroho et al. Reference Nugroho, Wang, Hosking and Milinazzo1999), where the moving load is prescribed. Here, the potential needs to be found and its solution requires the Green function. A closely related problem is that of a body advancing at forward speed
$U$
below a free surface and oscillating with frequency
$\omega$
. In the case of infinite water depth, when
$\tau =U\omega /g=1/4$
, where
$g$
is acceleration due to gravity, the corresponding Green function is singular and the solutions near the critical value were found to change very rapidly (Grue & Palm Reference Grue and Palm1985; Wu & Eatock Taylor Reference Wu and Eatock Taylor1987; Wu Reference Wu1991), and the problem at the critical value is not solvable. Later, Liu & Yue (Reference Liu and Yue1993) proved the solution for a real body of finite volume at
$\tau =1/4$
. They subsequently proposed a modified boundary integral equation (BIE) for the solution at this critical point. Palm & Grue (Reference Palm and Grue1999) further demonstrated that the solution for a foil of zero thickness with zero volume was also finite. However, the procedure of Liu & Yue (Reference Liu and Yue1993) relies on the condition that an introduced parameter
$\varGamma \neq 0$
, which may not always be straightforward to verify in more general cases. In this paper, we propose a novel procedure to solve the problem of a uniform current passing bodies submerged beneath an ice sheet at the critical points
$F=F_c$
and
$F=1$
. Although the detailed derivation and results are provided only for the 2-D case, the procedure can be used in 3-D problems. In particular, a new and modified BIE is derived in which the effect of the singularity is removed. This enables us to obtain the solution at the critical Froude number directly, which has not been achieved previously, and is highly significant. Moreover, the approach here is more general and straightforward, and it can be easily applied to a wide range of related problems.
The rest of the paper is arranged as follows. The boundary value problem and the Green function are introduced in § 2. In § 3, the procedure to treat singularity and a modified BIE are proposed. In § 4, case studies and results are presented for a single and double submerged circular cylinder, as well as an ellipse with different angles of attack, followed by the conclusions in § 5.
2. The boundary value problem and the Green function
The problem of a uniform current of speed
$U$
passing a body of arbitrary shape submerged below an ice sheet sketched in figure 1 is considered. A Cartesian coordinate system
$O\text{-}xz$
is defined, with the
$x$
-axis along the undisturbed water surface and against the direction of the current, and the
$z$
-axis pointing vertically upwards. The undisturbed water has depth
$H$
and its surface is covered by a homogeneous ice sheet.

Figure 1. Sketch of the problem and the defined coordinate system.
The fluid with density
$\rho$
is assumed to be inviscid, incompressible and homogeneous, and its motion is irrotational. The linearised velocity potential theory is employed for the problem, as in those works mentioned above. The total velocity potential is written as

where
$\phi$
denotes the disturbed velocity potential by the body, and is governed by the Laplace equation, or

The boundary condition on the ice sheet provides (Li et al. Reference Li, Wu and Shi2019)

where
$L$
represents the flexural rigidity of the ice sheet. On the body surface
$S_B$
and the seabed, the impermeable boundary condition gives


where
$\boldsymbol{n} = (n_x, n_z)$
represents the unit normal vector of
$S_B$
. Apart from that, the radiation condition should also be imposed at the far field, as a result of which the waves at
$x = +\infty$
and
$x = -\infty$
will have group velocities larger and smaller than
$U$
, respectively.
It is common that this kind of problem is solved through the BEM, in which the Green function
$G(\boldsymbol{P},\! \boldsymbol{Q})$
, or the velocity potential at the field point
$\boldsymbol{P}(x, z)$
due to a source located at
$\boldsymbol{Q}(x_0, z_0)$
is essential;
$G$
satisfies all the boundary conditions above, apart from that on the body surface
$S_B$
. From Li et al. (Reference Li, Wu and Shi2019), we have

where
$r_1 = \sqrt {(x - x_0)^2 + (z - z_0)^2}$
,
$r_2 = \sqrt {(x - x_0)^2 + (z + z_0 + 2H)^2}$
and



The integration route
$\mathcal{L}$
in (2.6) is from
$0$
to
$+\infty$
and its path at singularities corresponding to
$K(\alpha , U) = 0$
follows the procedure of Lighthill (Reference Lighthill1978). The nature of this equation has been discussed extensively by Yang, Wu & Ren (Reference Yang, Wu and Ren2024). When
$F_c \lt F = U / \sqrt {gH} \lt 1$
, where
$F_c$
denotes the critical Froude number (Davys et al. Reference Davys, Hosking and Sneyd1985),
$K(\alpha , U) = 0$
has two positive real roots
$\kappa _{-1}$
and
$\kappa _0$
with
$\kappa _{-1} \gt \kappa _0$
. The integration route
$\mathcal{L}$
should pass under (over) the first-order poles at
$\kappa _{-1}$
(
$\kappa _0$
). The integration then can be split into the Cauchy principal integration part and the residual part. Correspondingly there will be a
$\kappa _{-1}$
(
$\kappa _0$
) wave at
$x = +\infty$
(
$-\infty$
). When
$F \gt 1$
,
$\kappa _0$
becomes a purely negative imaginary root and its wave disappears. There is only one pole at
$\kappa _{-1}$
, and its wave at
$x = +\infty$
remains. When
$F \lt F_c$
,
$\kappa _{-1}$
and
$\kappa _0$
become a pair of conjugate complex roots, or
$\overline {\kappa }_{-1} = \kappa _0$
with
$\text{Re}\{\kappa _{-1}\} \gt 0$
and
$\text{Im}\{\kappa _{-1}\} \gt 0$
. There is no singularity on
$\mathcal{L}$
and no wave at infinity.
3. Solution procedure for submerged bodies at critical Froude numbers
With the help of the Green function, the differential equation for
$\phi$
can be converted into an integral equation over the body surface (Wehausen Reference Wehausen1973)

where
$\varLambda (\boldsymbol{P})$
denotes the 2-D solid angle at
$\boldsymbol{P}$
,
$\partial /\partial n_0$
represents the normal derivative at
$(x_0, z_0)$
. In general, the body surface can be discretised into small elements and (3.1) can be solved numerically. However, complexity arises at a critical speed. To demonstrate the problem explicitly, we may write

where
$R_c(\alpha ) \neq 0$
when
$\alpha \gt 0$
. Based on the discussion below (2.9), when
$F \to F_c$
,
$\kappa _0$
and
$\kappa _{-1}$
merge into
$\kappa _c$
and
$K(\alpha , U) \to (\alpha ^2 - \kappa _c^2)^2 R_c(\alpha )$
. In such a case, both
$K(\kappa _c, U_c) = K_\alpha (\kappa _c, U_c) = 0$
, where
$U_c = F_c \sqrt {gH}$
and
$K_\alpha$
denotes the derivative with respect to
$\alpha$
. As a result, the residue at
$\kappa _c$
becomes infinite, or the Green function at
$F_c$
is infinite, and (3.1) becomes unsolvable. This leads to the question of whether
$\phi$
is also infinite or exists. To answer this, we may define

where
$\mathcal{F}_c(\alpha , \boldsymbol{P},\! \boldsymbol{Q})$
does not contain any real pole. In such a case,
$G$
in (2.6) becomes

We may take out the singularities in (3.4), and write

Based on the definitions of
$\mathcal{L}$
,
$\kappa _0$
and
$\kappa _{-1}$
, using the residue theorem, we obtain

Substituting (3.5) and (3.6) into (3.4), we have

where


with
$\mathcal{F}_c (\alpha , \boldsymbol{P},\! \boldsymbol{Q})$
defined in (3.3). When
$F \to F_c$
,
$\left | \kappa _{-1} - \kappa _0 \right | \to 0$
, we may apply L’Hôpital’s rule to (3.8), which provides

Thus,
$\tilde {G}$
is finite at
$F = F_c$
, and is consistent with the result obtained by applying the Hadamard regularisation to the integral in (3.4). By contrast,
$\varPsi (\boldsymbol{P},\! \boldsymbol{Q})$
is singular. Besides, since
$(\kappa _{-1} - \kappa _0)$
is imaginary when
$F \lt F_c$
, but real when
$F_c \lt F \lt 1$
,
$\varPsi (\boldsymbol{P},\! \boldsymbol{Q})$
takes different forms when
$F \to F_c + 0^\pm$
.
We may substitute (3.7) into (3.1), which provides

The second integral in (3.11) is singular. To deal with this, we may investigate the nature of the Green function more closely. At large
$|x - x_0|$
, we may extend the integration route
$\mathcal{L}$
in (3.4) from
$\alpha \in (0, +\infty )$
to
$\alpha \in (-\infty , +\infty )$
. The path at the pole
$-\kappa _{-1}$
(
$-\kappa _0$
) is the same as that at
$\kappa _{-1}$
(
$\kappa _0$
) when it is real. When they are complex,
$-\kappa _{-1}$
(
$-\kappa _0$
) is in the lower (upper) half-plane. Additionally,
$\overline {\kappa }_{-1} = \kappa _0$
when
$F \lt F_c$
. In such a case, the residue theorem can be used in the upper or lower half-plane, depending on the sign of
$(x - x_0)$
. We need to keep only the leading
$\kappa _0$
and
$ $
$\kappa _{-1}$
terms, as well as the term due to the pole at
$\alpha = 0$
in
$G$
, because other poles are on the imaginary axis (Yang et al. Reference Yang, Wu and Ren2024) and decay far more rapidly. This provides


where the subscripts
$\pm \infty$
correspond to the sign of
$(x - x_0)$
,
$\mathscr{D}$
is a constant due to the pole at
$\alpha = 0$
and
$H(x)$
represents the Heaviside step function, reflecting the location of the pole at different Froude numbers. When
$F \lt 1$
, from (3.9), (3.12) and (3.13) we have

This gives

where
$\phi _{ \pm \infty }$
represents
$\phi$
at
$x \to \pm \infty$
. Applying (3.15) to (3.11), we have

Based on the behaviour of
$G_{\pm \infty }({\boldsymbol P},\! {\boldsymbol Q})$
at
$F \lt 1$
in (3.12) and (3.13), we may write


where
$A_\pm$
and
$B_\pm$
are unknowns,
$\mathscr{C}$
is known as the blockage constant for the free surface problem (Mei & Chen Reference Mei and Chen1976). When
$F \to F_c + 0^+$
, both
$\kappa _{-1}$
and
$\kappa _0$
are real and positive, and
$\phi _{\pm \infty }$
are wavy functions. By contrast, when
$F \to F_c + 0^-$
,
$\textrm{Im}\{\kappa _{-1}\} = -\textrm{Im}\{\kappa _0\} = \epsilon \to 0^+$
,
$\phi _{\pm \infty }$
decay extremely slowly with
$x$
. If
$F = F_c$
is taken,
$\epsilon = 0$
, and they become wavy functions too. Substituting (3.17) into (3.16), and letting
$F = F_c$
, we obtain

Compared with (3.11),
$\tilde {G}$
in (3.18) is finite. However, the equation has four additional unknowns, namely
$A_{\pm }$
and
$B_{\pm }$
. To resolve that, we may let
$|x|\to +\infty$
in (3.10), and invoke the theorem of residue for
$\tilde {G}$
. This provides

where
$\operatorname{sgn}(x-x_0)$
is the sign function,
$\mathscr{D}_c = g\pi /(gH - U_c^2)$
and

In (3.19),
$\tilde {G}_{\infty }({\boldsymbol P},\! {\boldsymbol Q})$
is associated with
$\kappa _c$
wave terms in (3.17), and
$\mathscr{D}_c |x - x_0|$
is linked to the blockage constant
$\mathscr{C}$
. Therefore, by letting
$|x| \to +\infty$
in (3.18), it becomes

Substituting (3.20) into (3.21), and matching both
$\sin (\kappa _c x)$
and
$\cos (\kappa _c x)$
terms of the equation, we obtain



where



Equation (3.22) then provides four additional equations. In such a case, the velocity potential
$\phi$
can be solved from (3.18) and (3.22). What is important here is that there is no singular term in these equations. Also, the same equations are valid for both
$F \to F_c + 0^\pm$
. Hence, the solution is continuous at
$F = F_c$
. This then resolves the difficulty caused by the singularity of the Green function at
$F = F_c$
.
Apart from the critical speed at
$F \to F_c$
, where
$(\kappa _{-1} - \kappa _0) \to 0$
, there is another critical speed at
$F \to 1$
, where
$\kappa _0 \to 0$
(Yang et al. Reference Yang, Wu and Ren2024). In particular, when
$F \to 1^-$
,
$\kappa _0 \to 0^+$
. As can be seen from (3.8) and (3.9), the Green function is singular. Using L’Hôpital’s rule in (3.3) at
$\alpha = \kappa _0$
, we have

Hence, (3.8) becomes

which is finite. Correspondingly (3.9) may be written as

where

Here,
$\varPsi (\boldsymbol{P},\! \boldsymbol{Q})$
is unbounded because of the
$1/\kappa _0$
term. The singular term may also be removed using the procedure in (3.14)–(3.16). In fact, at
$F = 1^-$
,
$\kappa _0 \to 0^+$
, we may apply the theorem of residue to (3.25). Letting
$|x| \to +\infty$
, we obtain

where
$\mathscr{D}_0 = ({9\pi }/{5H}) ( 1 - ( {5L}/{\rho g H^4}) )$
, and

We may further let
$F = 1^-$
in (3.17b
) and invoke L’Hôpital’s rule. This gives

Substituting (3.17a ) and (3.30) into (3.16), we have

We may perform an operation similar to that in (3.18). In (3.31), letting
$|x| \to +\infty$
, invoking (3.28) and matching terms of
$\sin (\kappa _{-1}x)$
,
$\cos (\kappa _{-1}x)$
,
$(z+H)^2 - x^2$
and
$x$
in the BIE. The following equations can be obtained:




where
$\mathcal{G}_j$
(
$j=1,2$
) are defined in (3.23a
) and (3.23b
), and



In such a case,
$A_\pm$
,
$B_\pm$
and
$\phi$
at
$F = 1^-$
can be solved from (3.31) and (3.32), and the solution is finite. Moreover, if we substitute (3.32) back into (3.31),
$\phi$
in fact can be solved directly from the BIE involving a modified Green function
$G^{\prime }(\boldsymbol{P},\! \boldsymbol{Q})$

where

This removes the singularity
$F=1^-$
and involves the unknowns only on the body surface. When
$F \to 1^+$
,
$\kappa _0 \to 0^- \times \textrm {i}$
(Yang et al. Reference Yang, Wu and Ren2024). We notice from (3.3) that
$\textrm{Re} \{ \mathcal{F}_c (\kappa _0, {\boldsymbol P},\! {\boldsymbol Q}) \} \to O(\kappa _0^{-1})$
, and thus (3.24) is no longer valid. To deal with that,
$\mathcal{F}_c (\kappa _0, {\boldsymbol P},\! {\boldsymbol Q})$
in (3.8) and (3.9) can be replaced with
$ [\mathcal{F}_c (\kappa _0, {\boldsymbol P},\! {\boldsymbol Q}) + \mathcal{F}_c (-\kappa _0, {\boldsymbol P},\! {\boldsymbol Q})]/2$
. Noticing

Hence,
$\tilde {G}({\boldsymbol P},\! {\boldsymbol Q})$
in (3.25) is still valid. The parameter
$\varPsi ({\boldsymbol P},\! {\boldsymbol Q})$
in (3.26) should be modified as

In such a case, (3.14)–(3.16) are still satisfied. However, since
$\kappa _0$
is imaginary here,
$G_{\pm \infty }$
in (3.12) and (3.13) contain the decaying wave component
$\kappa _0$
in both upstream and downstream regions. Based on the behaviour of
$G$
at infinity,
$\phi _{\pm \infty }$
should be rewritten as


where
$E_\pm$
are unknowns. Letting
$F=1^+$
in (3.38a
) and (3.38b
), we obtain


Following a similar procedure above, we may substitute (3.39) into (3.16). This provides

Letting
$|x| \to +\infty$
in (3.40), using (3.28), as well as matching functions
$\sin (\kappa _{-1} x)$
,
$\cos (\kappa _{-1} x)$
and
$x$
in the BIE, (3.32a
) and (3.32b
) can be also obtained for
$A_+$
and
$B_+$
. Besides, for
$E_\pm$
, we have


Hence, the velocity potential
$\phi$
at
$F = 1^+$
can be solved from (3.32a
), (3.32b
), (3.40) and (3.41). These equations have no singularity. The analysis above also shows that
$\phi$
at
$F = 1^\pm$
is discontinuous due to the difference in the modified equations.
4. Results and discussion
In the following computation, typical values of physical parameters of the ice sheet and the fluid domain are selected as (Li et al. Reference Li, Wu and Shi2019)

As shown in figure 2, two submerged body configurations are considered, namely, a single ellipse or two circular cylinders positioned beneath an ice sheet, for case studies. In the numerical implementation, each body surface is discretised into 128 linear elements, and equations are discretised and treated based on the procedure in Lu, He & Wu (Reference Lu, He and Wu2000), and the solution has been found to be convergent.
4.1. Verification through a single circular cylinder submerged below an ice sheet
Computation is first conducted for a single circular cylinder, which corresponds to a special case of figure 2(a) with dimensions
$b = a = H/8$
, and its centre is located at
$(x_c, z_c) = (0, -2a)$
. The distribution of the velocity potential
$\phi$
on the surface of the body is shown in figure 3, where
$\tan \theta = (z - z_c)/(x - x_c)$
. It can be seen that the results by the modified BIE at
$F = F_c$
and
$F = 1^{\pm }$
are fully consistent with those by the usual BIE, or (3.1), when
$F$
is sufficiently close to these critical Froude numbers. Besides, as shown in figures 3(b) and 3(c), noticeable differences can be observed in the curves of
$\phi$
versus
$\theta$
at
$F = 1^{-}$
and
$F = 1^{+}$
, which further verifies the proof in § 3, that
$\phi$
is discontinuous at
$F = 1$
. In such a case, the solution at
$F = 1$
is not unique, and any suitable linear combination of solutions at
$F \to 1^{-}$
and
$F \to 1^{+}$
can be a solution at
$F = 1$
. Furthermore, another interesting observation is that, when
$F \to 1^{+}$
, the solution of
$\phi$
tends to
$Ux$
. This can be confirmed by substituting it into (3.32a
), (3.32b
), (3.40) and (3.41), and it can be verified that all these equations are satisfied. It should be noted that
$\phi =Ux$
is the solution only for the case of
$F=1^+$
, as it does not satisfy the far-field conditions in any other cases, including
$F=1^-$
.

Figure 2. Sketch of the conducted case studies. (a) A single ellipse; (b) double circular cylinders.

Figure 3. Velocity potential around the surface of the single circular cylinder at critical speeds: (a)
$F\to F_c$
; (b)
$F\to 1^-$
; (c)
$F\to 1^+$
; (
$H=8a$
,
$(x_c,z_c )=(0,-2a)$
,
$F_c\approx 0.786882$
).

Figure 4. Resistance (a) and lift (b) on the circular cylinder versus the depth-based Froude number (
$H = 8a$
,
$(x_c, z_c) = (0, -2a)$
,
$F_c \approx 0.786882$
).
The resistance
$F_R$
and lift
$F_L$
on on the body can be evaluated by numerically integrating the pressure
$p$
over its surface, or

where

The minus sign in
$F_R$
means that the force is positive when it is in the direction of the incoming current. The corresponding results from BIE are shown in figure 4, where the solution based on the multipole expansion by Li et al. (Reference Li, Wu and Shi2019), applicable for
$F \neq F_c$
and
$F \neq 1$
, is also included for comparison and validation. Additionally, the hydrodynamic forces at
$F = F_c$
and
$F = 1^{\pm }$
are calculated using the modified BIE, with the numerical values being also presented in table 1. In figure 4(a),
$F_R = 0$
at
$F = F_c + 0^{-}$
, which can be known from the far-field formula (Yang et al. Reference Yang, Wu and Ren2021), since no wave exists at infinity. When
$F$
passes
$F = F_c$
, a rapid change occurs, followed by a peak in
$F_R$
at
$F \gt F_c$
. By contrast, in figure 4(b), a very large peak value of the lift
$F_L$
is observed before
$F = F_c$
, while at
$F = F_c$
,
$F_L$
is relatively small but non-zero. At
$F = 1^{\pm }$
, clear and sudden jumps in the hydrodynamic forces can be observed in both figures 4(a) and 4(b), indicating that the forces at this point are discontinuous. Besides, both
$F_R$
and
$F_L$
tend to zero as
$F \to 1^{+}$
since
$\phi \to Ux$
, as discussed above.
Table 1. Resistance and lift on the circular cylinder at critical Froude numbers (
$H = 8a$
,
$(x_c, z_c) = (0, -2a)$
,
$F_c \approx 0.786882$
).


Figure 5. Ice sheet deflection
$\eta$
due to a circular cylinder versus
$x$
as
$F \to F_c$
(
$H=8a$
,
$(x_c,z_c )=(0,-2a)$
,
$F_c\approx 0.786882$
).
The ice sheet deflection can be calculated from

Figure 5 shows the ice sheet deflection
$\eta (x)$
near the critical speed
$F = F_c$
. The maximum deflection occurs just before
$F = F_c$
, and forms a pronounced trough at
$F = 0.77$
. By contrast,
$\eta (x)$
is bounded and not particularly large at
$F = F_c$
. This is consistent with the variation trend of the lift observed in figure 4(b). This behaviour implies that subcritical speeds slightly smaller than
$F_c$
may be the most efficient speeds to generate sufficient force capable of fracturing the ice sheet. The parameter
$\eta (x)$
near
$F = 1^\pm$
is presented in figure 6. The results by the conventional BIE in (3.1) gradually tend to those by the modified BIE as
$F \to 1^\pm$
. However,
$\eta (x)$
tends to the result at
$F = 1^-$
much more quickly than that at
$F = 1^+$
. In particular, a significant difference can still be observed between results at
$F = 1.001$
and
$F = 1^+$
, even though the Froude numbers are very close. We may substitute (3.30) into (4.4). This provides
$\eta (x) \to U(-A_x + B_x)/g$
as
$F \to 1^-$
at
$x = -\infty$
, because the wavenumber downstream tends to
$0$
. By contrast, at
$F \to 1^+$
,
$\eta (x)$
becomes flat because
$\phi$
tends to
$Ux$
.

Figure 6. Ice sheet deflection
$\eta$
due to a circular cylinder versus
$x$
as
$F \to 1$
; (a)
$F\to 1^-$
; (b)
$F\to 1^+$
; (
$H=8a$
,
$(x_c,z_c )=(0,-2a)$
,
$F_c\approx 0.786882$
).
4.2. A uniform flow passing an ellipse submerged below an ice sheet
We next analyse a submerged elliptical cylinder with aspect ratio
$b/a = 2$
and its centre located at
$(x_c, z_c) = (0, -2a)$
, as shown in figure 2(a). Three angles of attack
$\beta = 0^\circ$
,
$15^\circ$
and
$30^\circ$
are considered as case studies for non-circular and asymmetric bodies. The lift
$F_L$
and resistance
$F_R$
against the depth-based Froude number
$F$
are shown in figure 7, which is obtained through the usual BIE. The values of
$F_L$
and
$F_R$
at
$F = F_c$
and
$F = 1^{\pm }$
are computed using the modified BIE and given in table 2. In figure 7, it can be seen that the variation trends of
$F_L$
and
$F_R$
are quite similar to those of a circular cylinder shown in figure 4. For
$F_R$
in figure 7(a), when
$F \gt 1$
, the effect of
$\beta$
on the resistance is significant. At a fixed value of
$F$
, as
$\beta$
increases, a clear increase in
$F_R$
can be seen. For
$F_L$
in figure 7(b), when
$F \lt F_c$
, at a fixed value of
$F$
,
$F_L$
increases with the angle of attack
$\beta$
. By contrast, when
$F \gt F_c$
, the influence of
$\beta$
on
$F_L$
is not so obvious. Additionally, the results in table 2 further indicate that
$F_L$
is bounded and relatively small at
$F = F_c$
, but a very significant peak occurs just before
$F = F_c$
. This again suggests that an effective strategy for ice breaking by an underwater vehicle may be to operate at speeds slightly below
$F = F_c$
, rather than exactly at
$F = F_c$
.
Table 2. Resistance and lift on an ellipse at critical Froude numbers (
$H = 8a$
,
$b = 2a$
,
$(x_c, z_c) = (0, -2a)$
,
$F_c \approx 0.786882$
).


Figure 7. Resistance (a) and lift (b) on the ellipse versus the depth-based Froude number under different angles of attack (
$H = 8a$
,
$b=2a$
,
$(x_c, z_c) = (0, -2a)$
,
$F_c \approx 0.786882$
).
4.3. A uniform flow passing double circular cylinders submerged below an ice sheet

Figure 8. Resistance (a) and lift (b) on the double circular cylinders versus the depth-based Froude number (
$H = 8a$
,
$z_1=z_2=-2a$
,
$d=4a$
,
$F_c \approx 0.786882$
).
A case study is also conducted for multiple submerged bodies, as shown in figure 2(b). We consider two identical circular cylinders with radius
$a=H/8$
, submerged at
$z_1=z_2=-2a$
, with the distance between their centres being
$d=4a$
(or
$x_1=-2a$
,
$x_2=2a$
). The indices 1 and 2 refer to the cylinders in the downstream and upstream regions, respectively. The values of
$F_L$
and
$F_R$
versus
$F$
are shown in figure 8, and the values at critical Froude numbers are shown in table 3. When
$F\lt F_c$
, similar to the case of a single circular cylinder (Li et al. Reference Li, Wu and Shi2019), the profile of the flexural gravity wave should be symmetrical about the origin, and there is no wave at the far field. Hence, the total resistance should be zero. However, as shown in figure 8(a), the resistance on each cylinder is non-zero. When
$F$
is small,
$F_R^{(1)}$
is positive and
$F_R^{(2)}$
is negative, which means mutual expulsion between the two bodies. As
$F$
increases and approaches
$F_c$
,
$F_R^{(1)}$
becomes negative, while
$F_R^{(2)}$
becomes positive, gradually creating a strong attraction effect between the two bodies. Different from a single body, it can be seen that the resistance has a large magnitude just below the critical speed. At
$F=F_c$
, the magnitudes of both resistances drop to relatively small values, which can be seen in table 3. As for lift in figure 8(b),
$F_L^{(1)}=F_L^{(2)}$
when
$F\lt F_c$
because of symmetry. When
$F_c\lt F\lt 1$
, both
$F_R^{(1)}$
and
$F_L^{(1)}$
are positive, whereas
$F_R^{(2)}$
and
$F_L^{(2)}$
are always negative. At
$F=1$
, clear jumps can also be found in the lifts and resistances in figure 8 and table 3. When
$F\gt 1$
, it can be seen that
$F_R^{(1)}$
remains positive, while
$F_R^{(2)}$
remains negative, and the attraction effect between two bodies becomes stronger as
$F$
increases within the range considered.
Table 3. Resistance and lift on two circular cylinders at critical Froude numbers (
$H = 8a$
,
$d = 4a$
,
$z_1 = z_2 = -2a$
,
$F_c \approx 0.786882$
).

5. Conclusion
The problem of a uniform current passing bodies submerged beneath an ice sheet at the critical Froude numbers
$F=F_c$
and
$F=1^\pm$
is investigated, based on the linearised velocity potential theory for the fluid and the elastic thin plate model for the ice sheet. It has been shown that, although the Green function is singular as
$F\to F_c$
and
$F\to 1$
, the velocity potential due to a real body remains finite. Particularly, the solution is continuous at
$F=F_c$
, bounded but discontinuous at
$F=1^\pm$
. Additionally, a modified BIE is proposed to solve
$\phi$
at these critical Froude numbers. The key to the success of this procedure has been to convert the singular terms in the conventional BIE to the velocity potential at the far field, thereby removing the singularity.
Various case studies are conducted to verify the mathematical proof and procedure. A single and double circular cylinder, as well as an elliptical cylinder with different angles of attack have been considered. The results are obtained through the usual BIE with
$F$
being sufficiently close to the critical Froude number and comparisons are made with results at the critical Froude number by the modified BIE. Detailed analyses are also conducted for the hydrodynamic forces on the submerged body and the generated flexural gravity waves. The following features have been observed. (i) The lift on the body is bounded and relatively small at
$F=F_c$
. However, when
$F$
is slightly below
$F_c$
, a large positive lift force occurs on the cylinder, which leads to a large deflection of the ice sheet. This may be used as an effective Froude number for ice breaking. (ii) When the depth-based Froude number
$F=1^+$
and
$F=1^-$
, clear differences are observed in both the forces and the ice sheet deflection profiles. This further indicates that the solution at
$F=1$
may not be unique. In particular, when
$F\to 1^+$
, it is found that
$\phi \to Ux$
, which leads the hydrodynamic forces on the body to tend to 0, and the deflection profile to become flat. (iii) For a single submerged body, the resistance is zero at
$F\lt F_c$
; For multiple submerged bodies, although the total resistance is zero at
$F\lt F_c$
, the resistance on each individual body is non-zero and it can be very large when the Froude number is just below the critical one. In the case of two submerged circular cylinders, depending on
$F$
, both attraction and repulsion effects can occur between the two bodies.
The present work has successfully resolved the challenge of solving the velocity potential problem through BIE with the Green function at the critical Froude numbers. This has allowed us to obtain some detailed results and gain some insights into the physics at and near the critical Froude number. What is more significant is that the solution procedure developed in this work is not just confined to the current problem, but can also be used in a wide range of related problems, for example, the singular problem in Liu & Yue (Reference Liu and Yue1993), the singular problem of fluid and structure interaction at critical Froude numbers or natural frequencies in an ice-covered channel Yang et al. (Reference Yang, Wu and Ren2021, Reference Yang, Wu and Ren2022) and also a uniform current interaction with a 3-D body in unbounded ocean.
Declaration of interests
The authors report no conflict of interest.