3.1. The Green function: velocity potential due to a single source

The Green function $G(x,y,z,x_0,y_0,z_0)$ is the velocity potential at point $P(x,y,z)$ due to a source at $P_0(x_0,y_0,z_0)$, which satisfies the following equation:

(3.1)\begin{equation} \nabla^2G=\delta(x-x_0)\delta(y-y_0)\delta(z-z_0), \end{equation}

where $\delta (x)$ is the Dirac delta function. Here $\xi (x,y,x_0,y_0,z_0)$ is defined as the wave elevation at point $(x,y)$ induced by the source at $P_0(x_0,y_0,z_0)$. Furthermore, $G(x,y,z,x_0,y_0,z_0)$ and $\xi (x,y,x_0,y_0,z_0)$ also need to satisfy the boundary conditions (2.6), (2.7) and (2.9)–(2.12).

Performing the Fourier transform for $G$ and $\xi$ in the $x$-direction

(3.2)\begin{equation} \left.\begin{gathered} \hat{G}=\int_{-\infty}^{+\infty}G\,\text{e}^{-\text{i}kx}\,{\text{d}\kern0.7pt x}\\ \hat{\xi}=\int_{-\infty}^{+\infty}\xi \,\text{e}^{-\text{i}kx}\,{\text{d}\kern0.7pt x} \end{gathered}\right\}, \end{equation}

and applying (3.2) to (3.1), we have

(3.3)\begin{equation} -k^2\hat{G}+\frac{\partial^2\hat{G}}{\partial y^2}+ \frac{\partial^2\hat{G}}{\partial z^2}=\delta(y-y_0)\delta(z-z_0)\,\text{e}^{-\text{i}kx_0}. \end{equation}

Based on the impermeable condition in (2.9), $\hat {G}$ can be further expanded into an orthogonal series of cosine functions in the $y$-direction. Using the condition in (2.10), the solution of (3.3) can be written in the following form:

(3.4)\begin{equation} \hat{G}=\sum_{n=0}^{+\infty}Z_n(k,z,x_0,y_0,z_0)\cos\sigma_n(y+b), \end{equation}

where

(3.5)\begin{align} & Z_n(k,z,x_0,y_0,z_0) \nonumber\\ &\quad ={-}\frac{[\exp({-K_n|z-z_0|})+\exp({-K_n(z+z_0+2H)})]\, \text{e}^{-\text{i}kx_0}\cos\sigma_n(y_0+b)}{2(1+\delta_{n0})bK_n} \nonumber\\ &\qquad +b_n\frac{\text{e}^{-\text{i}kx_0}\cosh K_n(z+H)}{\cosh K_nH}, \end{align}

$\sigma _n=n{\rm \pi} /2b$, $K_n=\sqrt {k^2+\sigma _n^2}$ and where $\delta _{ij}$ denotes the Kronecker delta function. The terms of $b_n$ in (3.5) correspond to the general solution of (3.3) when the right-hand side is zero. Here $b_n$ are to be determined by the conditions on the ice sheet. Applying the Fourier transform to (2.6) and (2.7), we have

(3.6)$$\begin{gather} (\rho g + Lk^4)\hat{\xi}-2Lk^2\frac{\partial^2 \hat{\xi}}{\partial y^2}+L \frac{\partial^4 \hat{\xi}}{\partial y^4}=\text{i}k\rho U\hat{G}, \end{gather}$$

(3.7)$$\begin{gather}-\text{i}kU\hat{\xi}=\frac{\partial \hat{G}}{\partial z}. \end{gather}$$

We choose to follow the approach taken by Ren et al. (Reference Ren, Wu and Li2020), which involves expanding $\partial ^4\hat {\xi }/\partial y^4$ rather than $\hat {\xi }$ into a cosine series, thus

(3.8)\begin{equation} \frac{\partial^4\hat{\xi}}{\partial y^4}=\text{e}^{-\text{i}kx_0}\sum_{n=0}^{+\infty}a_n\cos\sigma_n(y+b). \end{equation}

Then, through integration four times, $\hat {\xi }$ can be obtained as

(3.9)\begin{equation} \hat{\xi}=\text{e}^{-\text{i}kx_0}\left[c_0+c_1y+c_2y^2+c_3y^3+ \frac{a_0}{24}y^4+\sum_{n=1}^{+\infty}\frac{a_n}{\sigma_n^4}\cos\sigma_n(y+b)\right], \end{equation}

where $a_n\ (n=0,1,2\cdots )$ are unknown coefficients and are functions of $k$, $c_i\ (i=0\sim 3)$ are four constants which can be linked to $a_n$ through edge conditions. It should be mentioned here that $c_0$, $c_2$ and $a_{2n}$ correspond to symmetric components, while $c_1$, $c_3$ and $a_{2n+1}$ correspond to antisymmetric components. Substituting (3.4), (3.5) and (3.9) into (3.6) and (3.7), we have

(3.10)\begin{align} & (\rho g+Lk^4)\left[c_0+c_1y+c_2y^2+c_3y^3+\frac{a_0}{24}y^4+ \sum_{n=1}^{+\infty}\frac{a_n}{\sigma_n^4}\cos\sigma_n(y+b)\right] \nonumber\\ &\qquad -2k^2L\left[2c_2+6c_3y+\frac{a_0}{2}y^2-\sum_{n=1}^{+\infty} \frac{a_n}{\sigma_n^2}\cos\sigma_n(y+b)\right]+L\sum_{n=0}^{+\infty}a_n\cos\sigma_n(y+b) \nonumber\\ &\quad =\text{i}k\rho U\sum_{n=0}^{+\infty}\left[-\frac{\text{e}^{{-}K_nH}\cosh K_n(z_0+H)}{(1+\delta_{n0})bK_n}\cos \sigma_n(y_0+b)+b_n\right]\cos\sigma_n(y+b), \end{align}

(3.11)\begin{align} & -\text{i}kU\left[c_0+c_1y+c_2y^2+c_3y^3+\frac{a_0}{24}y^4+ \sum_{n=1}^{+\infty}\frac{a_n}{\sigma_n^4}\cos\sigma_n(y+b)\right] \nonumber\\ &\quad =\sum_{n=0}^{+\infty}\left[\frac{\text{e}^{{-}K_nH}\cosh K_n(z_0+H)}{(1+\delta_{n0})b} \cos\sigma_n(y_0+b)+b_nK_n\tanh K_nH\right]\cos\sigma_n(y+b). \end{align}

The term $y^j\ (\kern0.06em j=0\sim 4)$ can be further expanded into the orthogonal series of cosine functions as

(3.12)\begin{equation} y^j=\sum_{n=0}^{+\infty}d_n^{(\kern0.06em j)}\cos\sigma_n(y+b). \end{equation}

Then, (3.10) and (3.11) can be written as

(3.13)\begin{align} & (\rho g+Lk^4)\left[c_0d_n^{(0)}+c_1d_n^{(1)}+c_2d_n^{(2)}+c_3d_n^{(3)}+ \frac{a_0}{24}d_n^{(4)}+\frac{(1-\delta_{n0})a_n}{\sigma_n^4}\right] \nonumber\\ &\qquad -2k^2L\left[2c_2d_n^{(0)}+6c_3d_n^{(1)}+\frac{a_0}{2}d_n^{(2)}- \frac{(1-\delta_{n0})a_n}{\sigma_n^2}\right]+La_n \nonumber\\ &\quad =\text{i}k\rho U\left[-\frac{\text{e}^{{-}K_nH}\cosh K_n(z_0+H)}{(1+\delta_{n0})bK_n} \cos\sigma_n(y_0+b)+b_n\right],\quad n=0,1,2\cdots \end{align}

(3.14)\begin{align} & -\text{i}kU\left[c_0d_n^{(0)}+c_1d_n^{(1)}+c_2d_n^{(2)}+c_3d_n^{(3)}+ \frac{a_0}{24}d_n^{(4)}+\frac{(1-\delta_{n0})a_n}{\sigma_n^4}\right] \nonumber\\ &\quad =\left[\frac{\text{e}^{{-}K_nH}\cosh K_n(z_0+H)}{(1+\delta_{n0})b}\cos \sigma_n(y_0+b)+b_nK_n\tanh K_nH\right],\quad n=0,1,2\cdots . \end{align}

From (3.13) and (3.14), $a_n$ and $b_n$ can be expressed as

(3.15)$$\begin{gather} a_n = \alpha_{n,0}c_0+\alpha_{n,1}c_1+\alpha_{n,2}c_2+ \alpha_{n,3}c_3+\text{i}R_n,\quad n=0,1,2\cdots \end{gather}$$

(3.16)$$\begin{gather}b_n = \frac{1}{K_n\tanh K_nH}\left(\text{i}\beta_{n,0}c_0+\text{i}\beta_{n,1}c_1+\text{i}\beta_{n,2}c_2+\text{i}\beta_{n,3}c_3+S_n\right),\quad n=0,1,2\cdots \end{gather}$$

where

(3.17a)\begin{align} \alpha_{n,j} &= \frac{1}{\varDelta_n}\left\{\begin{aligned} & d_n^{(\kern0.06em j)}\left[(\rho g+Lk^4)K_n\tanh K_nH -\rho k^2U^2\right]\\ & -4\delta_{2j}k^2Ld_n^{(0)}K_n\tanh K_nH-12\delta_{3j}k^2Ld_n^{(1)}K_n\tanh K_nH \end{aligned}\right\} \nonumber\\ &\quad +(1-\delta_{n0})\gamma_n\alpha_{0,j}, \end{align}

(3.17b)\begin{equation} \beta_{n,j}={-}Uk\left[d_n^{(\kern0.06em j)}+\frac{d_n^{(4)}}{24}\alpha_{0,j}+ \frac{(1-\delta_{n0})}{\sigma_n^4}\alpha_{n,j}\right],\quad j=0\sim3 \end{equation}

and

(3.18a)$$\begin{gather} R_n=\frac{\rho Uk\cosh K_n(z_0+H)\cos\sigma_n(y_0+b)}{(1+\delta_{n0})b \varDelta_n\cosh K_nH}+(1-\delta_{n0})\gamma_nR_0, \end{gather}$$

(3.18b)$$\begin{gather}S_n=kU\left(\frac{d_n^{(4)}}{24}R_0+\frac{1-\delta_{n0}}{\sigma_n^4}R_n\right)-\frac{\text{e}^{{-}K_nH}\cosh K_n(z_0+H)\cos\sigma_n(y_0+b)}{(1+\delta_{n0})b}, \end{gather}$$

with

(3.19a)$$\begin{gather} \varDelta_n=\left\{\begin{aligned} & \delta_{n0}\left[-\left(\frac{\rho g+Lk^4}{24}d_n^{(4)}-k^2Ld_n^{(2)}+L\right)K_n\tanh K_nH+ \frac{\rho k^2U^2d_n^{(4)}}{24}\right]\\ & \quad +(1-\delta_{n0})\left[-\left(\frac{\rho g+Lk^4}{\sigma_n^4}+ \frac{2k^2L}{\sigma_n^2}+L\right)K_n\tanh K_nH+\frac{\rho k^2U^2}{\sigma_n^4}\right] \end{aligned}\right\}, \end{gather}$$

(3.19b)$$\begin{gather}\gamma_n=\frac{1}{\varDelta_n}\left[\left(\frac{\rho g+Lk^4}{24}d_n^{(4)}-k^2Ld_n^{(2)}\right)K_n \tanh K_nH-\frac{\rho k^2U^2d_n^{(4)}}{24}\right]. \end{gather}$$

If we consider the symmetric or antisymmetric nature of $y^j$ $(\kern0.06em j=0\sim 3)$, we will have $d_{2n}^{(2j+1)}=d_{2n+1}^{(2j)}=0$ $(\kern0.06em j=0, 1)$. This leads to $\gamma _{2n+1}=0$ and $\alpha _{2n+1,2j}=\alpha _{2n,2j+1}=0$ $(\kern0.06em j=0, 1,\ n\geq 0)$. As a result, $a_n$ can be further expressed as

(3.20)\begin{equation} \left.\begin{gathered} a_{2n}=\alpha_{2n,0}c_0+\alpha_{2n,2}c_2+\text{i}R_{2n}\\ a_{2n+1}=\alpha_{2n+1,1}c_1+\alpha_{2n+1,3}c_3+\text{i}R_{2n+1} \end{gathered}\right\}. \end{equation}

It can be seen from (3.20) that $a_{2n}$ depend on two unknown coefficients $c_0$ and $c_2$, while $a_{2n+1}$ depend on $c_1$ and $c_3$; $b_n$ can be also treated in a similar way. The four unknown coefficients $c_i$ $(i=0\sim 3)$ can be determined from the four edge conditions, including those in (2.11). As a result, a system of linear equations of the following form can be established:

(3.21)\begin{equation} [\boldsymbol{\mathcal{A}}][\boldsymbol{C}]=[\boldsymbol{\mathcal{B}}], \end{equation}

where $[\boldsymbol {\mathcal {A}}]$ is a $4\times 4$ coefficient matrix, $[\boldsymbol {C}]$ is a column containing $c_i$, and $[\boldsymbol {\mathcal {B}}]$ is a known column. For a specific case, this $4\times 4$ matrix equation may be further simplified. If the edge conditions at $y=\pm b$ are same, (3.21) may be split into two independent $2\times 2$ submatrix equations, one for $c_0$ and $c_2$, and the other for $c_1$ and $c_3$. As a result, the symmetric and antisymmetric transverse waves in (3.9) become independent to each other. In general, $c_i$ $(i=0\sim 3)$ are fully coupled to each other, which leads $a_{2n}$ and $a_{2n+1}$ in (3.20) to become coupled. The solution of $c_i$ can be written as

(3.22)\begin{equation} c_i ={-}\text{i}\sum_{m=0}^{+\infty}\frac{\rho Uk\cosh K_m(z_0+H)\cos \sigma_m(y_0+b)}{(1+\delta_{m0})b|\boldsymbol{\mathcal{A}}|\varDelta_m\cosh K_mH}c_{m,i}^{\prime},\quad i=0\sim3, \end{equation}

where $|\boldsymbol {\mathcal {A}}|$ is the determinant of $[\boldsymbol {\mathcal {A}}]$. The elements of $[\boldsymbol {\mathcal {A}}]$ and $[\boldsymbol {\mathcal {B}}]$, as well as coefficient $c_{m,i}^\prime$ are related to edge conditions, an example of these under clamped–clamped edges are given in Appendix A.

Once $c_i$ $(i=0\sim 3)$, $a_n$ and $b_n$ $(n=0,1,2\cdots )$ are obtained, then $G$ and $\xi$ can be obtained through the inverse Fourier transform. Using (Linton & McIver Reference Linton and McIver2001)

(3.23a) \begin{align} & \int_{-\infty}^{+\infty}\frac{\delta_{n0}\exp(-K_nH)-\exp({-K_n|z-z_0|+\text{i}k(x-x_0)})}{2K_n}\text{d} k \nonumber\\ &\quad =\delta_{n0}\ln \left(\frac{r}{H}\right)-(1-\delta_{n0})\mathscr{K}_0(\sigma_nr), \end{align}

(3.23b) \begin{align} & \int_{-\infty}^{+\infty}\frac{\delta_{n0}\exp(-K_nH)-\exp({-K_n(z+z_0+2H)+\text{i}k(x-x_0)})}{2K_n}\text{d} k \nonumber\\ &\quad =\delta_{n0}\ln \left(\frac{r^{\prime}}{H}\right)-(1-\delta_{n0})\mathscr{K}_0(\sigma_nr^\prime), \end{align}

where $r=\sqrt {(x-x_0)^2+(z-z_0)^2}$ and $r^{\prime }=\sqrt {(x-x_0)^2+(z+z_0+2H)^2}$, $\mathscr {K}_n$ denotes the $n$th-order modified Bessel function of second kind. We have the Green function

(3.24)$$\begin{gather} G(x,y,z,x_0,y_0,z_0)=\sum_{n=0}^{+\infty}\mathcal{Z}_n(x,z,x_0,y_0,z_0)\cos\sigma_n(y+b), \end{gather}$$

(3.25)$$\begin{gather}\mathcal{Z}_n(x,z,x_0,y_0,z_0)=\mathcal{Z}_{n,1}+\mathcal{Z}_{n,2}, \end{gather}$$

and

(3.26a) \begin{align} \mathcal{Z}_{n,1} &= \frac{1}{4{\rm \pi} b}\left\{\delta_{n0}\left[\ln \left(\frac{r}{H}\right)+\ln\left(\frac{r^{\prime}}{H}\right)\right]-2(1-\delta_{n0}) \left[\mathscr{K}_0(\sigma_nr)+ \mathscr{K}_0(\sigma_n r^{\prime})\right]\right\} \nonumber\\ &\quad \times\cos\sigma_n(y_0+b), \end{align}

(3.26b) \begin{align} \mathcal{Z}_{n,2}& =\frac{\rho U}{2{\rm \pi} b}\int_{\mathscr{L}}\sum_{m=0}^{+\infty}\frac{\begin{array}{l}\left\{ k\mu_{n,m}\cosh K_m(z_0+H) \cos\sigma_m(y_0+b)\right.\\\left. \quad \times\left[\cosh K_n(z+H)\exp({\text{i}k(x-x_0)})- \delta_{n0}\right]\right\}\end{array}}{(1+\delta_{m0})|\boldsymbol{\mathcal{A}}|\varDelta_mK_n\sinh K_nH\cosh K_mH}\text{d} k, \end{align}

where

(3.27)\begin{align} \mu_{n,m} &= \sum_{j=0}^{3}\beta_{n,j}c_{m,j}^{\prime}+\delta_{m0}|\boldsymbol{\mathcal{A}}|Uk \left(\frac{d_n^{(4)}}{24}+\frac{(1-\delta_{n0})\gamma_n}{\sigma_n^4}\right) \nonumber\\ &\quad +\delta_{nm}|\boldsymbol{\mathcal{A}}|\left[\frac{Uk(1-\delta_{n0})}{\sigma_n^4}- \frac{\varDelta_n(1+\text{e}^{{-}2K_nH})}{2\rho Uk}\right]. \end{align}

It should be mentioned that a constant term is, respectively, added to $\mathcal {Z}_{n,1}$ and $\mathcal {Z}_{n,2}$ in (3.26) to remove the high-order singularity at $k=0$, which will not affect the results as all the equations for $G$ involve only its spatial derivatives. The term $\mathcal {Z}_{n,1}$ is obtained based on the derivation in Li, Wu & Ren (Reference Li, Wu and Ren2021). There will be singularities in $\mathcal {Z}_{n,2}$ when $|\boldsymbol {\mathcal {A}}|(k)=0$. Because $|\boldsymbol {\mathcal {A}}|(k)$ is an even function, all its real roots can be represented by $\pm k_s$ $(s=1\cdots S)$, with $k_s>0$. In the case of Ren et al. (Reference Ren, Wu and Li2020) for wave propagation in the channel, each $k_s$ corresponds to the dispersion relationship between the wave frequency and wavenumber. Here, mathematically, each $k_s$ corresponds to a singularity in the integrand of the inverse Fourier transform. Physically, it corresponds to each wave generated by the source in the channel. The number of singularities can be more than one, and the value of $S$ depends on the current speed $U$ when other parameters are fixed, which will be discussed in detail later. The way to treat each singularity will depend on the group velocity of its corresponding wave. The wave will be at upstream $(x=+\infty )$ or downstream $(x=-\infty )$ if its group velocity is larger or smaller than $U$, respectively. Based on this, the integration route $\mathscr {L}$ in $\mathcal {Z}_{n,2}$ from $-\infty$ to $+\infty$ can be defined as follows. We may consider the integration of $f(k)\,\text {e}^{\text {i}kx}/(k-k_s)$ along the path $\mathscr {L}$. This can be split into the principle integration and a contribution from the pole. If the group velocity of the wave component $k_s$ is larger (smaller) than $U$, $\mathscr {L}$ passes over (under) the singularities at ${\pm }k_s$. Thus, the contribution from the pole at $k=k_s$ will be $-\text {i}{\rm \pi} f(k_s)\,\text {e}^{\text {i}k_sx}$ or $+\text {i}{\rm \pi} f(k_s)\,\text {e}^{\text {i}k_sx}$. If we use (Wehausen & Laitone Reference Wehausen and Laitone1960)

(3.28)\begin{equation} \lim_{|x|\to+\infty}\text{P.V.}\int_{-\infty}^{+\infty} \frac{f(k)}{k-k_s}\text{e}^{\text{i}kx}\,\text{d} k=\text{sgn}(x)\text{i}{\rm \pi} f(k_s)\,\text{e}^{\text{i}k_sx}, \end{equation}

in the integral in $\mathcal {Z}_{n,2}$, we can find that the radiation condition (2.12) is satisfied.

It may be of interest to see that the Green function is symmetric about the source and field points, or $G(x,y,z,x_0,y_0,z_0)=G(x_0,y_0,z_0,x,y,z)$, which is shown in Appendix B.

3.2. Multipole expansion for the horizontal circular cylinder

The Green function derived above can be used to convert the governing equation to an integral equation over the surface of a body with arbitrary shape. For some special geometries, such as a sphere, the solution may be found through expansion in terms of the multipole obtained through differentiating the Green function with respect to the position of the source (see, for example, Wu (Reference Wu1998b)). For a horizontal circular cylinder, the potential can be expanded into the cosine series as used for $G$. The governing Laplace equation in $(x,y,z)$ then becomes the modified Helmholtz equation in $(x,z)$ for a circular section. Subsequently, the solution of the modified Helmholtz equation can be obtained from a two-dimensional multipole expansion. To construct that, we may apply a source distribution $\varsigma (y_0)$ along the centre line of the cylinder. This is then expanded into the cosine series, or $\varsigma (y_0)=\sum _{n=0}^{+\infty }\mathcal {V}_n\cos \sigma _n(y_0+b)$. It will create the following potential:

(3.29)\begin{equation} \varphi=\sum_{n=0}^{+\infty}{\mathcal{V}_n\int_{{-}b}^{b}G\cos{\sigma_n(y_0+b)}\,\text{d} y_0}, \end{equation}

where $G$ is the Green function derived in the previous section. Substituting (3.24) into (3.29), we obtain

(3.30)\begin{equation} \varphi=\sum_{n=0}^{+\infty}{\varphi_n\cos{\sigma_n\left(y+b\right)}}, \end{equation}

where

(3.31)\begin{equation} \varphi_n=\varphi_{n,1}+\varphi_{n,2}, \end{equation}

with

(3.32a)\begin{align} \varphi_{n,1}&=\frac{\mathcal{V}_n}{2{\rm \pi}}\left\{\delta_{n0} \left[\ln \left(\frac{r}{H}\right)+\ln\left(\frac{r^{\prime}}{H}\right)\right]-(1-\delta_{n0}) \left[\mathscr{K}_0(\sigma_n r)+\mathscr{K}_0(\sigma_n r^{\prime})\right]\right\}, \end{align}

(3.32b)\begin{align} \varphi_{n,2}&=\frac{\rho U}{2{\rm \pi}} \nonumber\\ &\quad \times\int_{\mathscr{L}}\sum_{n^\prime=0}^{+\infty}\frac{\mathcal{V}_{n^\prime} \mu_{n,n^\prime}k\cosh K_{n^\prime}(z_0+H)\left[\cosh K_n(z+H)\exp({\text{i}k(x-x_0)})-\delta_{n0}\right]}{|\boldsymbol{\mathcal{A}}|\varDelta_{n^\prime}K_n \sinh K_nH \cosh K_{n^\prime}H}\text{d} k. \end{align}

We may use the operator

(3.33)\begin{equation} \left(D_{{\pm}}\right)^m={-}\frac{1}{2^{m-1}\left(m-1\right)!} \left(\frac{\partial}{\partial z_0}\pm \text{i}\frac{\partial}{\partial x_0}\right)^m. \end{equation}

This gives (Linton & McIver Reference Linton and McIver2001)

(3.34a)$$\begin{gather} \left(D_{{\pm}}\right)^m\ln{r}=\frac{\text{e}^{{\pm} \text{i}m\theta}}{r^m}, \end{gather}$$

(3.34b)$$\begin{gather}\left(D_{{\pm}}\right)^m\mathscr{K}_0\left(\sigma_nr\right)= \frac{\left({-}1\right)^m\sigma_n^m}{2^{m-1}\left(m-1\right)!}\text{e}^{{\pm} \text{i}m\theta}\mathscr{K}_m\left(\sigma_nr\right), \end{gather}$$

where $z-z_0=r\cos {\theta }$ and $x-x_0=r\sin {\theta }$. As mentioned in Li, Wu & Shi (Reference Li, Wu and Shi2019), because $\varphi$ is a real function, we may apply only the operator $(D_+)^m$. Using

(3.35a)$$\begin{gather} \left(D_+\right)^m\exp({\pm K_nz_0\pm \text{i}kx_0})={-}\frac{({\pm} 1)^m(K_n-k)^m}{2^{m-1}(m-1)!}\exp({\pm K_nz_0\pm \text{i}kx_0}), \end{gather}$$

(3.35b)$$\begin{gather}\left(D_+\right)^m\exp({\pm K_nz_0\mp \text{i}kx_0})={-}\frac{({\pm} 1)^m(K_n+k)^m}{2^{m-1}(m-1)!}\exp({\pm K_nz_0\mp \text{i}kx_0}), \end{gather}$$

as well as (3.34), the velocity potential of the multipole can be expressed as

(3.36)$$\begin{gather} \left(\varphi_+\right)^m=\left(D_{+}\right)^m\varphi= \sum_{n=0}^{+\infty}\left(\varphi_+\right)_n^m\cos{\sigma_n(y+b)}, \end{gather}$$

(3.37)$$\begin{gather}\left(\varphi_+\right)_n^m=\left(\varphi_+\right)_{n,1}^m+\left(\varphi_+\right)_{n,2}^m+ \left(\varphi_+\right)_{n,3}^m, \end{gather}$$

where $(\varphi _+)_{n,1}^m+(\varphi _+)_{n,2}^m=(D_{+})^m\varphi _{n,1}$, $(\varphi _+)_{n,3}^m=(D_{+})^m\varphi _{n,2}$ and can be obtained as

(3.38a)$$\begin{gather} \left(\varphi_+\right)_{n,1}^m=\frac{\mathcal{V}_n}{2{\rm \pi}} \left[\delta_{n0}\frac{\text{e}^{\text{i}m\theta}}{r^m}-(1-\delta_{n0}) \frac{({-}1)^m\sigma_n^m}{2^{m-1}(m-1)!}\text{e}^{\text{i}m\theta}\mathscr{K}_m(\sigma_nr)\right], \end{gather}$$

(3.38b)$$\begin{gather}\left(\varphi_+\right)_{n,2}^m=\frac{({-}1)^m\mathcal{V}_n}{2^{m+1}(m-1)!{\rm \pi}} \int_{-\infty}^{+\infty}\frac{(K_n-k)^m\exp({-K_n(z+z_0+2H)+\text{i}k(x-x_0)})}{K_n}\text{d} k, \end{gather}$$

(3.38c)\begin{gather} \hspace{-21pc}(\varphi_+)_{n,3}^{m} = \frac{-\rho U}{2^{m+1}(m-1)!{\rm \pi}} \nonumber\\ \hspace{3.5pc}\times\int_{\mathscr{L}}\sum_{n^\prime=0}^{+\infty} \frac{\mathcal{V}_{n^\prime}\mu_{n,n^\prime}k\cosh K_n(z+H)E_{n^{\prime},m}(k,z_0)\exp({\text{i}k(x-x_0)})}{|\boldsymbol{\mathcal{A}}|\varDelta_{n^\prime}K_n \sinh K_nH \cosh K_{n^{\prime}}H}\text{d} k, \end{gather}

with

(3.39)\begin{align} E_{n^{\prime},m}(k,z_0)=(K_{n^{\prime}}+k)^m\exp({K_{n^{\prime}}(z_0+H)})+ ({-}1)^m(K_{n^{\prime}}-k)^m\exp({-K_{n^{\prime}}(z_0+H)}). \end{align}

The potential $\phi$ due to the cylinder then can be written in a multipole expansion form as

(3.40)\begin{equation} \phi=\text{Re}\left\{\sum_{m=1}^{+\infty}{r_0^mg_m(\varphi_+)^m}\right\}. \end{equation}

It then satisfies all the boundary conditions met by the Green function and the only remaining one is that on the cylinder surface. To satisfy the body surface boundary condition, we may write the potential in the polar coordinate system. Using (Abramowitz & Stegun Reference Abramowitz and Stegun1970)

(3.41)\begin{equation} \left.\begin{gathered} \exp({K_n(z-z_0)\pm \text{i}k(x-x_0)})=\sum_{m=0}^{+\infty}T_{n,m}(r) \left[A_{n,m}^\pm(k)\,\text{e}^{\text{i}m\theta}+A_{n,m}^\mp(k)\,\text{e}^{-\text{i}m\theta}\right]\\ \exp({-K_n(z-z_0)\pm \text{i}k(x-x_0)})=\sum_{m=0}^{+\infty}({-}1)^mT_{n,m}(r) \left[A_{n,m}^\mp(k)\,\text{e}^{\text{i}m\theta}+A_{n,m}^\pm(k)\,\text{e}^{-\text{i}m\theta}\right] \end{gathered}\right\}, \end{equation}

where

(3.42)\begin{equation} \left.\begin{gathered} A_{n,m}^+(k)=\delta_{n0}k^m+\left(1-\delta_{n0}\right)\left(\frac{K_n+k}{\sigma_n}\right)^m\\ A_{n,m}^-(k)=\left(1-\delta_{n0}\right)\left(1-\delta_{m0}\right)\left(\frac{\sigma_n}{K_n+k}\right)^m\\ T_{n,m}(r)=\delta_{n0}\frac{r^m}{m!}+\left(1-\delta_{n0}\right)\mathscr{I}_m\left(\sigma_nr\right) \end{gathered}\right\}, \end{equation}

and $\mathscr {I}_m$ denotes the $m$th-order modified Bessel function of first kind, we have

(3.43) \begin{equation} \phi=\text{Re}\left\{\sum_{n=0}^{+\infty}\sum_{m=0}^{+\infty} \begin{aligned} & \left[\begin{aligned} & Q_{n,m}f_{n,m}\,\text{e}^{\text{i}m\theta}+T_{n,m} \sum_{m^\prime=1}^{+\infty}f_{n,m^\prime}\left(\mathcal{C}_{n,m,m^\prime}^+ \text{e}^{\text{i}m\theta}+\mathcal{C}_{n,m,m^\prime}^-\text{e}^{-\text{i}m\theta}\right)\\ & \quad +T_{n,m}\sum_{n^\prime=0}^{+\infty}\sum_{m^\prime=1}^{+\infty}f_{n^\prime,m^\prime} \left(\mathcal{D}_{n,n^\prime,m,m^\prime}^+\text{e}^{\text{i}m\theta}+ \mathcal{D}_{n,n^\prime,m,m^\prime}^-\text{e}^{-\text{i}m\theta}\right) \end{aligned}\right]\\ &\times\cos{\sigma_n(y+b)} \end{aligned}\right\}, \end{equation}

where

(3.44a)$$\begin{gather} Q_{n,m}(r)=\frac{(1-\delta_{m0})}{2{\rm \pi}}\left[\delta_{n0}\left(\frac{r_0}{r}\right)^m- \frac{(1-\delta_{n0})(-\sigma_nr_0)^m}{2^{m-1}(m-1)!}\mathscr{K}_m(\sigma_nr)\right], \end{gather}$$

(3.44b)$$\begin{gather}\mathcal{C}_{n,m,m^{\prime}}^{{\pm}}=\frac{({-}1)^{m+m^{\prime}} r_0^{m^{\prime}}}{2^{m^{\prime}+1}(m^{\prime}-1)!{\rm \pi}}\int_{-\infty}^{+\infty} \frac{(K_n-k)^{m^{\prime}}\exp({-2K_n(z_0+H)})A_{n,m}^{{\mp}}}{K_n}\text{d} k, \end{gather}$$

(3.44c)\begin{gather} \hspace{-19pc}\mathcal{D}_{n,n^{\prime},m,m^{\prime}}^{{\pm}}=\frac{-\rho Ur_0^{m^{\prime}}}{2^{m^{\prime}+2}(m^{\prime}-1)!{\rm \pi}} \nonumber\\ \hspace{2pc} \times\int_{\mathscr{L}} \frac{\mu_{n,n^\prime}k\left[\exp({K_n(z_0+H)})A_{n,m}^{{\pm}}+({-}1)^m\exp({-K_n(z_0+H)})A_{n,m}^{{\mp}}\right] E_{n^{\prime},m^{\prime}}}{\varDelta_{n^\prime}|\boldsymbol{\mathcal{A}}|K_n \sinh K_nH \cosh K_{n^{\prime}}H}\text{d} k.\nonumber\\ \end{gather}

It should be mentioned here that a new unknown coefficient $f_{n,m}$ is defined as $f_{n,m}=\mathcal {V}_ng_m$. The impermeable condition on the body surface in (2.8) gives

(3.45)\begin{align} & Q_{n,m}^\prime(r_0)f_{n,m}+T_{n,m}^\prime(r_0) \sum_{m^\prime=1}^{+\infty}\left(f_{n,m^\prime}\mathcal{C}_{n,m,m^\prime}^+{+}{\bar{f}}_{n,m^\prime}{\bar{\mathcal{C}}}_{n,m,m^\prime}^-\right) \nonumber\\ &\quad +T_{n,m}^\prime(r_0)\sum_{n^\prime=0}^{+\infty}\sum_{m^\prime=1}^{+\infty} \left(f_{n^\prime,m^\prime}\mathcal{D}_{n,n^\prime,n,m^\prime}^+{+}{\bar{f}}_{n^\prime,m^\prime}{\bar{\mathcal{D}}}_{n,n^\prime,m,m^\prime}^-\right)={-}\text{i} \delta_{n0}\delta_{m1}U, \end{align}

where $n=0,1,2\cdots$ and $m=1,2,3\cdots$ $Q_{n,m}^\prime$ ($T_{n,m}^\prime$) represents the derivative of $Q_{n,m}$ ($T_{n,m}$) with respect to $r$, which can be obtained as

(3.46)\begin{equation} \left.\begin{gathered} Q_{n,m}^\prime(r)=\frac{(1-\delta_{m0})}{2{\rm \pi}}\left[-\delta_{n0}m\frac{r_0^m}{r^{m+1}}+ \frac{(1-\delta_{n0})(-\sigma_nr_0)^{m+1}}{2^{m-1}(m-1)!r_0}\mathscr{K}_m^\prime(\sigma_nr)\right]\\ T_{n,m}^\prime(r)=\delta_{n0}\frac{r^{m-1}}{(m-1)!}+(1-\delta_{n0})\sigma_n\mathscr{I}_m^\prime(\sigma_nr) \end{gathered}\right\}. \end{equation}

After (3.45) is solved, $\phi$ can then be obtained.

3.3. Hydrodynamic forces

Once the potential is found, the lift $F_L$ and resistance $F_R$ on the cylinder can be calculated through the integration of hydrodynamic pressure over the body surface. Thus, we have

(3.47)\begin{equation} -\text{i}F_R+F_L={-}\iint_{S_B}p\,\text{e}^{\text{i}\theta}\,\text{d} S. \end{equation}

The pressure $p$ can be obtained by the Bernoulli equation

(3.48)\begin{equation} p={-}\tfrac{1}{2}\rho\boldsymbol{\nabla}(\phi-Ux)\boldsymbol{\cdot}\boldsymbol{\nabla}(\phi-Ux). \end{equation}

We may notice that the product term is kept here, as it may be small on the ice sheet but may not be on the body surface. When determining the gradient term in (3.48) in the cylindrical coordinate system, we may substitute (3.45) into (3.43) and have

(3.49)\begin{equation} \left.\begin{gathered} \left.\frac{\partial(\phi-Ux)}{\partial r}\right|_{r=r_0}=0\\ \left.\frac{\partial(\phi-Ux)}{\partial\theta}\right|_{r=r_0}=\text{Re}\left\{\sum_{n=0}^{+\infty} \sum_{m=0}^{+\infty}{\text{i}m\psi_{n,m}\,\text{e}^{\text{i}m\theta}\cos\sigma_n(y+b)}\right\}\\ \left.\frac{\partial(\phi-Ux)}{\partial y}\right|_{r=r_0}={-}\text{Re}\left\{\sum_{n=0}^{+\infty} \sum_{m=0}^{+\infty}{\sigma_n\psi_{n,m}\,\text{e}^{\text{i}m\theta}\sin\sigma_n(y+b)}\right\} \end{gathered}\right\}, \end{equation}

where

(3.50)\begin{equation} \psi_{n,m}=f_{n,m}\left[Q_{n,m}(r_0)-Q_{n,m}^\prime(r_0)\frac{T_{n,m}(r_0)}{T_{n,m}^\prime(r_0)}\right]. \end{equation}

Substituting (3.48) and (3.49) into (3.47), we obtain

(3.51)\begin{equation} -\text{i}F_R+F_L=\frac{{\rm \pi}\rho b}{2r_0}\sum_{n=0}^{+\infty} \sum_{m=0}^{+\infty}\left[(1+\delta_{n0})m(m+1)+(1-\delta_{n0}) \sigma_n^2r_0^2\right]\psi_{n,m}{\bar{\psi}}_{n,m+1}. \end{equation}

Similar to that in Wu (Reference Wu1995), the resistance can be also obtained by a far-field integral, which is shown in Appendix C.

3.4. Deflection of the ice sheet

We may use (2.7) to obtain the expression of $\eta$, and this gives

(3.52)\begin{equation} \eta(x,y)={-}\frac{1}{U}\int{\frac{\partial\phi}{\partial z}{\text{d}\kern0.7pt x}}+C(y). \end{equation}

Substituting (3.36)–(3.40) into (3.52), we have

(3.53)$$\begin{gather} \eta=\text{Im}\left\{\sum_{m=1}^{+\infty}\sum_{n=0}^{+\infty}{r_0^m\eta_{n,m}} \cos\sigma_n(y+b)\right\}+C(y), \end{gather}$$

(3.54)$$\begin{gather}\eta_{n,m}=\eta_{n,m}^{(1)}+\eta_{n,m}^{(2)}+\eta_{n,m}^{(3)}, \end{gather}$$

where

(3.55a)$$\begin{gather} \eta_{n,m}^{(1)}=\frac{f_{n,m}}{2^{m+1}(m-1)!{\rm \pi} U}\int_{-\infty}^{+\infty}\frac{(K_n+k)^m\exp({K_nz_0+\text{i}k(x-x_0)})}{k}\text{d} k, \end{gather}$$

(3.55b)$$\begin{gather}\eta_{n,m}^{(1)}=\frac{({-}1)^mf_{n,m}}{2^{m+1}(m-1)!{\rm \pi} U}\int_{-\infty}^{+\infty}\frac{(K_n-k)^m\exp({-K_n(z_0+2H)+\text{i}k(x-x_0)})}{k}\text{d} k, \end{gather}$$

(3.55c)$$\begin{gather}\eta_{n,m}^{(3)}=\frac{\rho}{2^{m+1}(m-1)!{\rm \pi}}\int_{\mathscr{L}} \sum_{n^{\prime}=0}^{+\infty}\frac{f_{n^{\prime},m}\mu_{n,n^\prime}E_{n^{\prime},m}(k,z_0) \exp({\text{i}k(x-x_0)})}{\varDelta_{n^{\prime}}|\boldsymbol{\mathcal{A}}|\cosh K_{n^{\prime}}H}\text{d} k \end{gather}$$

and where $C(y)$ in (3.52) is the integration constant. As $C(y)$ is not function of $x$, we may determine it at $x\to +\infty$. Based on (3.36)–(3.40) and (3.53)–(3.55), the asymptotic expressions of $\phi$ and $\eta$ at $x\to +\infty$ can be written as

(3.56)$$\begin{gather} \phi=\text{Re}\left\{\sum_{j=1}^{S}\phi^{(\kern0.06em j)}(y,z)\exp({-\text{i}k_jx})\right\}+\text{sgn}(x)\phi^{(0)}(y,z), \end{gather}$$

(3.57)$$\begin{gather}\eta=\frac{1}{U}\text{Re}\left\{\sum_{j=1}^{S}\frac{1}{\text{i}k_j}\frac{\partial \phi^{(\kern0.06em j)}(y,0)}{\partial z}\exp({-\text{i}k_jx})\right\}+C(y), \end{gather}$$

where $\phi ^{(0)}$ is due to the singularity of the integrand at $k=0$ and is related to the ‘blockage parameter’ (Mei & Chen Reference Mei and Chen1976). It should be noted that the summation in (3.56) and (3.57) should include only those terms with group velocity larger than $U$. It can be shown that $\partial \phi ^{(0)}/\partial z=0$ on $z=0$ and therefore it does not contribute to $\eta$. We may combine (2.6) and (2.7) to eliminate $\eta$, and then use (3.56) in the result. We have

(3.58)\begin{align} & \left[L\left(k_j^4-2k_j^2\frac{\partial^2}{\partial y^2}+ \frac{\partial^4}{\partial y^4}\right)+\rho g\right] \frac{\partial\phi^{(\kern0.06em j)}}{\partial z}-\rho U^2k_j^2\phi^{(\kern0.06em j)}=0, \nonumber\\ & \quad z=0,\quad x\to+\infty,\enspace j=1\sim S. \end{align}

Substituting (3.56) and (3.57) into (2.6), and using (3.58), we obtain

(3.59)\begin{equation} L\frac{\text{d}^4C}{{\text{d}y}^4}+\rho gC=0. \end{equation}

The boundary conditions of $C$ can be established by substituting (3.57) into (2.11) which is satisfied by the first term on the right-hand side of (3.57). This gives, at $y=\pm b$,

(3.60)\begin{equation} \left.\begin{gathered} C=0,\quad C_y=0,\quad \text{clamped edge}\\ C=0,\quad C_{yy}=0,\quad \text{simply supported edge}\\ C_{yy}=0,\quad C_{yyy}=0,\quad \text{free edge} \end{gathered}\right\}. \end{equation}

It can be found that $C(y)=0$ under any of these edge conditions.

3.5. Numerical procedure

Although the expression for the potential has been written explicitly, some of the computations still have to be performed numerically. Taking into account that the integrand decays exponentially, the integration for $k$ from $-\infty$ to $+\infty$ is truncated at a sufficiently large value $K_T$ and treated as

(3.61)\begin{align} \text{P.V.}\int_{-\infty}^{+\infty}F(k)\,\text{d} k=\text{P.V.}\int_{0}^{+\infty} \left[F(k)+F({-}k)\right]\text{d} k\approx\text{P.V.}\int_{0}^{K_T}\left[F(k)+F({-}k)\right]\text{d} k, \end{align}

where P.V. indicates the Cauchy principal value. The range $(0,K_T)$ is divided into many small steps and the Gaussian method is used in each step for integration. To deal with multiple singularities in the integrand, the following numerical procedures are applied. Assume that $F(k)$ contains $n$ first-order singularities in $k\in (0,+\infty )$. We may write

(3.62)\begin{equation} F(k)=\frac{G(k)}{\prod_{i=1}^{n}(k-k_i)}=G(k)\sum_{i=1}^{n} \frac{1}{\prod_{j=1(\kern0.06em j\neq i)}^{n}(k_i-k_j)}\times\frac{1}{k-k_i}, \end{equation}

its Cauchy principal value can be calculated through

(3.63)\begin{align} \text{P.V.}\int_{0}^{+\infty}F(k)\,\text{d} k&=\sum_{i=1}^{n}\frac{1}{\prod_{j=1(\kern0.06em j\neq i)}^{n}(k_i-k_j)} \nonumber\\ &\quad \times\left[\int_{0}^{2k_i}\frac{G(k)-G(k_i)}{k-k_i}\text{d} k+ \int_{2k_i}^{+\infty}\frac{G(k)}{k-k_i}\text{d} k\right], \end{align}

where $G(k)$ can be calculated by

(3.64)\begin{equation} \left.\begin{gathered} G(k)=F(k)\prod_{i=1}^{n}(k-k_i),\quad k\neq k_i\\ G(k_i)=\prod_{j=1(\kern0.06em j\neq i)}^{n}(k_i-k_j)\lim_{k\to k_i}(k-k_i)F(k),\quad i=1\sim n \end{gathered}\right\}. \end{equation}

By doing that, all the singularity effects are eliminated. Although the original one integral is split into $n$ integrals, it is still more computationally efficient, as the integration step ${\rm \Delta} k$ can be appropriately chosen for each integral. For a single integral of $F(k)$, ${\rm \Delta} k$ has to be very small, especially when some of the singularities are very close to each other.

When determining the residue at the singularities caused by $|\boldsymbol {\mathcal {A}}|(k_s)=0$, the value of $|\boldsymbol {\mathcal {A}}|^\prime (k_s)$ needs to be calculated numerically. Here we adopt the method proposed by Li et al. (Reference Li, Wu and Ren2021). Assuming $f(k)/(k-k_s)=g(k)/|\boldsymbol {\mathcal {A}}|(k)=P(k)$, we have

(3.65)\begin{equation} f(k_s)=\frac{g(k_s)}{|\boldsymbol{\mathcal{A}}|^\prime(k_s)}= \lim_{k\to k_s}(k-k_s)P(k)\approx\frac{P(k_s+{\rm \Delta} k)-P(k_s-{\rm \Delta} k)}{2}{\rm \Delta} k. \end{equation}

4.1. Analysis of the dispersion relationship of an ice-covered channel

As discussed after (3.24), there are singularities in the integrand when $|\boldsymbol {\mathcal {A}}|(k,U)=|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$ (where $\hat {k}=kH$), which are equivalent to the dispersion relationship and correspond to the waves at infinity. Compared with the two-dimensional case of a homogeneous ice plate with infinite extent (Li et al. Reference Li, Wu and Shi2019), the dispersion relationship here is more complicated. The root $Fn$ of the equation at each given $\hat {k}$ is not unique. In fact, there is an infinite number of solutions of $Fn$ at a given $\hat {k}$, we denote these solutions as $Fn^{(i)}$ $(i=1,2,3\cdots )$, with $Fn^{(1)}< Fn^{(2)}< Fn^{(3)}<\cdots$. Similar to that in Khabakhpasheva et al. (Reference Khabakhpasheva, Shishmarev and Korobkin2019), the curves $Fn^{(i)}$ against $\hat {k}$ are shown in figure 2. It can be seen that for each curve there is a minimum value at $\hat {k}={\hat {k}}_c^{(i)}$, which can be called $i$th-order critical Froude number and be denoted by $Fn_c^{(i)}$. The value of each $Fn_c^{(i)}$ can be found in figure 2. When $Fn< Fn_c^{(1)}$, there will be no solution in $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$, or there will be no waves propagating to $x=\pm \infty$. When $Fn_c^{(1)}< Fn< Fn_c^{(2)}$, there will be two solutions in $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$. The ones corresponding to $\hat {k}<{\hat {k}}_c^{(1)}$ and $\hat {k}>{\hat {k}}_c^{(1)}$ will lead to waves at $x=-\infty$ and $x=+\infty$, respectively, because their group velocities are smaller and larger than $U$, respectively, as in Li et al. (Reference Li, Wu and Shi2019). In fact, the non-dimensionalized group velocity $c_g^{(s)}$ can be obtained as

(4.2)\begin{equation} c_g^{\left(s\right)}=Fn({\hat{k}}_s)+{\hat{k}}_s\frac{\text{d}}{\text{d}\hat{k}}Fn({\hat{k}}_s), \end{equation}

where $({\hat {k}}_s, Fn({\hat {k}}_s))$ is an intersection point of $Fn$ with curve $\hat {k}-Fn^{(i)}$. It can be observed from figure 2 and (4.2) that when ${\hat {k}}_s<{\hat {k}}_c^{(i)}$ $({\hat {k}}_s>{\hat {k}}_c^{(i)})$, the slope $\text {d} Fn(\hat {k}_s)/\text {d}\hat {k}<0$ $(\text {d} Fn(\hat {k}_s)/\text {d}\hat {k}>0)$, and the group velocity $c_g^{(s)}< Fn(\hat {k}_s)$ $(c_g^{(s)}>Fn(\hat {k}_s))$ and wave will be at $x=-\infty$ $(x=+\infty )$.

Figure 2. Purely positive solution of $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$ for the clamped–clamped edge condition.

As $Fn$ increases, it will reach the second critical point $Fn_c^{(2)}$. Beyond that, there will be two more solutions in $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$. This leads to two more waves and one each at $x=-\infty$ and $x=+\infty$. In general, when $Fn_c^{(i)}< Fn< Fn_c^{(i+1)}$, there will be $2i$ roots in $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$, and $i$ waves at $x=-\infty$ and $x=+\infty$, respectively. However, the curve $Fn_c^{(1)}$ is different from the others. When $\hat {k}\to 0$, we can see that $\lim _{\hat {k}\to 0}Fn^{(1)}=Fn_0$ while $\lim _{\hat {k}\to 0}Fn^{(i)}=+\infty$ $(i\geq 2)$. Therefore, when $Fn>Fn_0$, corresponding to curve $Fn^{(1)}$, there will be only one solution in $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)=0$, and there will be $2i-1$ solutions when $Fn_c^{(i)}< Fn< Fn_c^{(i+1)}$.

As mentioned in § 3.1, the symmetric and antisymmetric transverse waves are completely independent in the cases of symmetric edges. If we further write $|\boldsymbol {\mathcal {A}}|$ as in (A10), it can be found that all the intersection points on curve $Fn^{(2i-1)}$ $(i\geq 1)$ in figure 2 are solutions of $|\boldsymbol {\mathcal {A}^S}|(\hat {k},Fn)=0$ and correspond to a symmetric mode, while those on curve $Fn^{(2i)}$ are solutions of $|\boldsymbol {\mathcal {A}^A}|(\hat {k},Fn)=0$ and correspond to an antisymmetric mode.

In the two-dimensional case, there is only one critical $Fn_c$, below which there will be no wave, and above which will be two waves at $x=-\infty$ and $x=+\infty$, respectively. When $Fn>1$, the wave at the downstream region will disappear but the one at the upstream region will remain. Here, when there is an infinite number of critical Froude numbers, $Fn_c^{(1)}$ is similar to $Fn_c$ in the two-dimensional case, below which there will be no wave. As $Fn$ increases and passes through each $Fn^{(i)}$ $(i=1,2,3\cdots )$ each pair of waves will be created, with one at the downstream and the other at the upstream regions. Here, $Fn_0$ here corresponds to $Fn=1$ in the two-dimensional case. When $Fn>Fn_0$, the wave at $x=-\infty$ due to $Fn^{(1)}$ will disappear, which is similar to the two-dimensional case. However, different from the two-dimensional case, there will be still waves at the downstream region.

The feature of wave in infinite water depth, or $H\to +\infty$, is also investigated, other parameters in (4.1) are kept the same. For the two-dimensional case, the critical value $Fn=1$ will disappear if $H\to +\infty$. There will always be two waves, with one each at the downstream and upstream regions, respectively, when the current speed is larger than the critical value. Here, using $\lim _{H\to +\infty } \tanh K_nH = 1$, $\alpha _{n,j}$ in (3.17a) can be further expressed as

(4.3)\begin{equation} \alpha_{n,j}=\frac{1}{\varDelta_n^{\infty}}\left\{\begin{aligned} & d_n^{(\kern0.06em j)}[(\rho g+Lk^4)K_n-\rho k^2U^2]\\ & \quad -4k^2Ld_n^{(0)}K_n\delta_{2j}-12k^2Ld_n^{(1)}K_n\delta_{3j} \end{aligned}\right\}+(1-\delta_{n0})\gamma_n^\infty\alpha_{0j}, \end{equation}

where

(4.4a)\begin{gather} \varDelta_n^{\infty}=\left\{\begin{aligned} & \delta_{n0}\left[-\left(\frac{\rho g+Lk^4}{24}d_n^{(4)}-k^2Ld_n^{(2)}+L\right)K_n+ \frac{\rho k^2U^2}{24}d_n^{(4)}\right]\\ & \quad +(1-\delta_{n0})\left[-\left(\frac{\rho g+Lk^4}{\sigma_n^4}+ \frac{2k^2L}{\sigma_n^2}+L\right)K_n+\frac{\rho k^2U^2}{\sigma_n^4}\right] \end{aligned}\right\}, \end{gather}

(4.4b)\begin{gather} \gamma_n^{\infty}=\frac{1}{\varDelta_n^\infty} \left[\left(\frac{\rho g+Lk^4}{24}d_n^{\left(4\right)}-k^2Ld_n^{\left(2\right)}\right)K_n- \frac{\rho k^2U^2}{24}d_n^{\left(4\right)}\right]. \end{gather}

If the half-width $b$ is chosen as the characteristic length, two new parameters can be defined as $\widetilde {Fn}=U/\sqrt {gb}$ and $\tilde {k}=kb$. Substituting (4.3) into (A2), the solution of $|\boldsymbol {\mathcal {A}}|(k,U)=|\boldsymbol {\mathcal {A}}|(\tilde {k},\widetilde {Fn})=0$ at $H\to +\infty$ can be obtained and are shown in figure 3. Similar to $Fn_c^{(i)}$ in figure 2, the $i$th-order critical value at $\tilde {k}={\tilde {k}}_c^{(i)}$ is defined as ${\widetilde {Fn}}_c^{(i)}$. It can be seen that the curve ${\widetilde {Fn}}^{(1)}$ is different from $Fn^{(1)}$ in figure 2. When $\tilde {k}\to 0$, we can find that $\lim _{\tilde {k}\to 0}\widetilde {Fn}^{(1)}=+\infty$. Therefore, there will always be $2i$ solutions in $|\boldsymbol {\mathcal {A}}|(\tilde {k},\widetilde {Fn})=0$ when ${\widetilde {Fn}}_c^{(i)}<\widetilde {Fn}<{\widetilde {Fn}}_c^{(i+1)}$, which leads to $i$ waves at $x=-\infty$ and $x=+\infty$, respectively.

Figure 3. Purely positive solution of $|\boldsymbol {\mathcal {A}}|(\tilde {k},\widetilde {Fn})=0$ for clamped–clamped edges at $H\to +\infty$.

We may also investigate how the above results vary with the channel width $b$, ice sheet thickness $h_i$ and various edge constraints while other values in (4.1) remain unchanged. Here $Fn_c^{(i)}$ $(i=1,2,3,4)$ at different $b$ and $h_i$ under clamped–clamped edges is shown in figure 4. It can be seen from figure 4(a) that all the $Fn_c^{(i)}$ decrease when $b$ increases, and they all tend to the value of $Fn_c$ without the tank wall, or the two-dimensional case (Li et al. Reference Li, Wu and Shi2019). In other words, when $b$ increases, the curves $Fn^{(i)}$ in figure 2 will gradually move towards each other and approach the curve in the two-dimensional case. Figure 4(b) shows that at $h_i=0$ all the $Fn_c^{(i)}=0$, which is expected as this is a free surface problem. As $h_i$ increases, all the $Fn_c^{(i)}$ increase. As $h_i$ becomes very large or $L$ is very large, the ice sheet becomes a rigid plate and there will be no waves, which means $Fn_c^{(i)}\to +\infty$. In figures 5(a) and 5(b), $Fn_c^{(1)}$ under other edge conditions show variation trends similar to that in figure 4, apart from the free–free one where $Fn_c^{(1)}$ increases slightly with $b$. Plots of $Fn_0$ versus with $b$ and $h_i$ are given in figures 5(c) and 5(d), respectively. It can be shown that under the free–free condition, $Fn_0=1$, which can be confirmed by substituting $\hat {k}=0$ and $Fn=1$ in $|\boldsymbol {\mathcal {A}}|(\hat {k},Fn)$. This is the same as the two-dimensional case. Under other edge conditions $Fn_0$ is different from $1$. It varies with $b$ in a way similar to that of $Fn_c^{(1)}$. It can be observed that at sufficiently large $b$, $Fn_0\to 1$, which is consistent with the two-dimensional case. Also $\lim _{h_i\to 0}{Fn_0}=1$, which is consistent with the free surface problem. Another interesting feature in figure 5 is when all the parameters are fixed, $Fn_c^{(1)}$ and $Fn_0$ increase in the following sequence: free–free, simply supported–free, clamped–free, simply supported–simply supported, clamped–simply supported and clamped–clamped.

Figure 4. Critical Froude number under the clamped–clamped edge condition: (a) $Fn_c^{(i)}$ $(i=1\sim 4)$ at different channel widths; (b) $Fn_c^{(i)}$ $(i=1\sim 4)$ at different ice sheet thickness.

Figure 5. Here $Fn_c^{(1)}$ and $Fn_0$ under different edge conditions: (a) $Fn_c^{(1)}$ at different channel widths; (b) $Fn_c^{(1)}$ at different ice thickness; (c) $Fn_0$ at different channel widths; (d) $Fn_0$ at different ice thickness. The clamped edge is denoted by C, the simply supported edge is denoted by SS and the free edge is denoted by F.

4.2. Deflection of the ice sheet due to an underwater source

Numerical results are first given for the ice sheet deflection $\xi$ induced by a submerged source. The quiescent water depth $H$ is used as the characteristic length. Here, $\xi$ can be obtained directly by using inverse Fourier transformation

(4.5)\begin{equation} \xi = \frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}\hat{\xi}\,\text{e}^{\text{i}kx}\,\text{d} k. \end{equation}

Substituting (3.7) into (4.5) and using (3.24)–(3.27), we have

(4.6)\begin{align} & \xi = \frac{\text{i}\rho}{2{\rm \pi} b}\nonumber\\ &\quad \times\int_{\mathscr{L}}\sum_{n=0}^{+\infty} \sum_{m=0}^{+\infty}\frac{\lambda_{n,m}\cosh K_m(z_0+H)\exp({\text{i}k(x-x_0)}) \cos\sigma_n(y+b)\cos\sigma_m(y_0+b)}{(1+\delta_{m0})\varDelta_m|\boldsymbol{\mathcal{A}}|\cosh K_mH}\text{d} k, \end{align}

where

(4.7)\begin{equation} \lambda_{n,m}=\sum_{j=0}^{3}\beta_{n,j}c_{m,j}^{\prime}+ \delta_{m0}|\boldsymbol{\mathcal{A}}|Uk\left(\frac{d_n^{(4)}}{24}+ \frac{(1-\delta_{n0})\gamma_n}{\sigma_n^4}\right)+\delta_{nm}(1-\delta_{n0}) \frac{Uk|\boldsymbol{\mathcal{A}}|}{\sigma_n^4}. \end{equation}

It should be noted that the infinite series in (4.6) is truncated at $n=m=15$ in the calculation, and the results have been found to have converged.

The ice deflection along the $y$-axis induced by a source at two different locations under clamped–clamped edges is shown in figure 6. In figure 6(a), the source is located at $(x_0,y_0,z_0)=(0,0,-H/5)$. This corresponds to a symmetric case about $y=0$, or only waves corresponding to even $m$ in (4.6) will exist. As $x=0$ and $x_0=0$, the principal integral of $\xi$ in (4.6) is zero. The remaining waves are those components which will propagate to $x=+\infty$ and $x=-\infty$, respectively. It can be seen from figure 6(a) that the wave profiles at $Fn=0.9$ and $Fn=1.05$ are dominated by the first symmetric transverse mode ($\sigma _2$), and at $Fn=1.4$ and $Fn=1.8$ the second symmetric mode ($\sigma _4$) become important, further, at $Fn=2.4$, the third symmetric mode ($\sigma _6$) also becomes important. In general, when $Fn_c^{(2m-1)}< Fn< Fn_c^{(2m+1)}$ $(m\geq 1)$, the contribution of the $m$th symmetric mode $(\sigma _{2m})$ to the overall wave profile will be significant. As $c_1=c_3=0$ and $a_{2n-1}=b_{2n-1}=0$ in (3.5) and (3.9), the singularities corresponding to $|\boldsymbol {\mathcal {A}^A}|(\hat {k},Fn)=0$ (see (A 6)) in the solution for these coefficients are no longer relevant. Therefore, the curves of $Fn^{(2i)}(\hat {k})$ can be removed from figure 2 in such a case. The transverse wave profiles with the source located at $(x_0,y_0,z_0)=(0,-b/2,-H/5)$ are shown in figure 6(b), which is an asymmetric case. In such a case, all the longitudinal wave components of $k_j$ $(\kern0.06em j=1\sim S)$ in (4.6) exist, and all the transverse waves of $\sigma _m$ will also exist. The overall transverse wave in this case is asymmetric.

Figure 6. Ice sheet deflection $\xi$ along the $y$-axis due to a source when the ice sheet is under the clamped–clamped edge condition: (a) $(x_0,y_0,z_0)=(0,0,-H/5)$; (b) $(x_0,y_0,z_0)=(0,-b/2,-H/5)$. The Froude numbers in the figure are in the following order: $Fn_c^{(1)}=0.8180<0.9< Fn_c^{(2)}=0.9981<1.05< Fn_c^{(3)}=1.3177<1.4< Fn_c^{(4)}=1.7682<1.8< Fn_c^{(5)}=2.3230<2.4< Fn_c^{(6)}=2.9600$.

The hydroelastic waves along the $x$-axis due to a source at $(x_0,y_0,z_0)=(0,0,-H/5)$ are given in figure 7. It can be seen from figure 7(a) that there is no wave propagating to $x=\pm \infty$ when $Fn< Fn_c^{(1)}$. The longitudinal wave profile is antisymmetric about $x=x_0$, which can be confirmed by the sine function in the principal integral of (4.6). As $Fn$ increases, it can be seen that in figures 7(a), 7(c) and 7(d) that the wave elevation near $x=x_0$ is significantly increased when $Fn$ is near each of $Fn_c^{(i)}$ $(i=1,3,5\cdots )$. When $Fn\to Fn_0$, the smaller root of $Fn^{(1)}(\hat {k})=Fn$ tends to zero, and the corresponding wavelength at $x=-\infty$ tends to infinity, which can be observed from figure 7(b). Since waves corresponding to the roots of $Fn^{(2i)}=Fn$ will not appear in the symmetric case, there is no wave propagating to $x=-\infty$ when $Fn_0< Fn< Fn_c^{(3)}$. As discussed in Appendix B, due to the net mass outflow from the source, or, mathematically, the singularity of function $\hat {\xi }$ at $k=0$, opposite mean wave elevations can be observed at $x\to +\infty$ and $x\to -\infty$, which are denoted by ${\bar {\xi }}_{+\infty }$ and ${\bar {\xi }}_{-\infty }$, respectively. It is also interesting to see that ${\bar {\xi }}_{+\infty }>0$ (${\bar {\xi }}_{-\infty }<0$) when $Fn< Fn_0$, while ${\bar {\xi }}_{+\infty }<0$ (${\bar {\xi }}_{-\infty }>0$) when $Fn>Fn_0$. In figure 8, the ice sheet deflection along $y=-b/2$ generated by a source located at $(x_0,y_0,z_0)=(0,-b/2,-H/5)$ is plotted. Phenomena similar to those in figure 7 may be also observed, apart from that waves corresponding to the roots of $Fn^{(2i)}=Fn$ also appear. In fact, from figures 8(c) and 8(d), we can see clearly the wave behind the source due to $Fn^{(2)}$ when $Fn>Fn_c^{(2)}$.

Figure 7. Ice sheet deflection $\xi$ along $x$-axis due to a source located at $(x_0,y_0,z_0)=(0,0,-H/5)$ when the ice sheet is under the clamped–clamped edge condition: (a) $Fn< Fn_c^{(1)}$; (b) $Fn_c^{(1)}< Fn< Fn_0$; (c) $Fn_c^{(2)}< Fn< Fn_c^{(3)}$; (d) $Fn_c^{(4)}< Fn< Fn_c^{(5)}$; where $Fn_c^{(1)}=0.8180$, $Fn_c^{(2)}=0.9981$, $Fn_c^{(3)}=1.3177$, $Fn_c^{(4)}=1.7682$, $Fn_c^{(5)}=2.3230$, $Fn_0=1.1469$.

Figure 8. Ice sheet deflection $\xi$ along $y=-b/2$ due to a source located at $(x_0,y_0,z_0)=(0,-b/2,-H/5)$ when the ice sheet is under the clamped–clamped edge condition: (a) $Fn< Fn_c^{(1)}$; (b) $Fn_c^{(1)}< Fn< Fn_c^{(2)}$; (c) $Fn_c^{(2)}< Fn< Fn_0$; (d) $Fn_0< Fn< Fn_c^{(3)}$; where $Fn_c^{(1)}=0.8180$, $Fn_c^{(2)}=0.9981$, $Fn_c^{(3)}=1.3177$, $Fn_0=1.1469$.

4.3. Hydrodynamic forces on a submerged horizontal cylinder

A submerged horizontal cylinder with its centre line located at $x_0=0$ and $z_0=-H/5$ in the current is considered next. The radius of the cylinder $r_0$ is employed as the characteristic length. For the free surface problem, the two-dimensional steady solution may no longer be possible within some ranges of the Froude number at a given submergence (Haussling & Coleman Reference Haussling and Coleman1979; Scullen & Tuck Reference Scullen and Tuck1995; Semenov & Wu Reference Semenov and Wu2020). This is mainly because the steady wave amplitude reaches its limit and may break, which is reflected by the fact the velocity magnitude at the wave peak tends to zero (Semenov & Wu Reference Semenov and Wu2020). This can occur when the submergence of the circular cylinder is chosen as $|z_0 |/r_0=2$ for the free surface flow. However, here the problem is quite different. On the one hand, the covering ice sheet will obviously affect the behaviour of the waves. More importantly, the case considered here is a three-dimensional one. For clamped edges as an example, the wave elevation along the edge is in fact zero. Indeed, it also has been found from the obtained numerical results that the magnitude of the velocity on the far-field is always positive for the case considered, apart from that near the critical Froude numbers. This suggests that the linear model is valid apart from when the Froude number is near one of the critical values. At the critical value, the results of the linear model become inaccurate, while its prediction for the critical value itself may be accurate. In such a case, nonlinearity and other effects may need to be taken into account.

To conduct numerical computations, the infinite series in (3.45) is truncated at ${n^\prime =N}$ and $m^\prime =M$. A convergence test is then carried out for $N$ and $M$ through the hydrodynamic forces on the cylinder when the ice sheet is under the clamped–clamped edge condition. The far-field formula (C7) is also employed to check the resistance $F_R$ obtained by the near field formula (3.51). The convergence test results are shown in figure 9. An excellent agreement can be seen between the results obtained by ${N=15}$ and $M=15$, and $N=25$ and $M=25$, which means the convergence has been achieved. Therefore, $N=15$ and $M=15$ is applied in the following calculation. We can also observe that the resistance calculated by (C7) is graphically the same as that from (3.51). It can be seen from figure 9(a) that $F_R=0$ when $Fn< Fn_c^{(1)}$, which is similar to the two-dimensional case considered by Li et al. (Reference Li, Wu and Shi2019). As mentioned above, there is no wave at infinity when $Fn< Fn_c^{(1)}$, which leads to the fact that coefficients $Q_{n,m}$, $\mathcal {C}_{n,m,m^\prime }^\pm$ and $\mathcal {D}_{n,n^\prime ,m,m^\prime }^\pm$ in (3.45) are all real, and correspondingly the solution $f_{n,m}$ should be imaginary. Together with (3.50), $\psi _{n,m}$ is imaginary, then from (3.51), we have $F_R=0$.

Figure 9. Convergence study for the hydrodynamic forces on the cylinder when the ice sheet is under the clamped–clamped edge condition: (a) resistance; (b) lift.

We then consider the resistance and lift under three different symmetric edge constraints. The results are provided in figure 10. In figure 10(a), when the Froude number passes $Fn=Fn_c^{(i)}$ $(i\geq 1)$ or $Fn_0$, the resistance changes rapidly. The number of the wave components at far-field will also change. In particular the resistance $F_R$ is zero initially. When $Fn$ passes $Fn=Fn_c^{(1)}$, $F_R$ jumps up rapidly. As $Fn$ passes $Fn_0$, $F_R$ will have a sharp drop. These are similar to the behaviours of the resistance at $Fn=Fn_c$ and $Fn=1$ in the two-dimensional case considered by Li et al. (Reference Li, Wu and Shi2019). When $Fn$ continues to increase and pass $Fn_c^{(2i-1)}$ $(i\geq 2)$, $F_R$ will first increase rapidly and then decrease rapidly, and return to almost the original trajectory of the curve. In the two-dimensional case, critical values of $Fn_c^{(i)}$ $(i>1)$ do not exist. For the $F_L$ given in figure 10(b), it is non-zero even at small $Fn$ because the flow is not symmetric about $z=z_0$. The lift near $Fn_c^{(2i-1)}$ $(i\geq 1)$ and $Fn_0$ have a behaviour similar to that of resistance. At large $Fn$ the lift may become negative similar to the free surface problem (Wu & Eatock Taylor Reference Wu and Eatock Taylor1987). It is also interesting to see from figure 10 that the differences in $F_R$ or $F_L$ under various edge conditions are much more obvious in the region near $Fn_c^{(i)}$ and $Fn_0$. When $Fn$ passes $Fn_c^{(2i-1)}$ $(i>1)$, we can see that the jump phenomena of $F_L$ and $F_R$ in the case of the free–free edges are not as obvious as those of the clamped–clamped and simply supported–simply supported edges. The variation of $F_R$ and $F_L$ against $Fn$ under three asymmetric edges constraints are shown in figure 11. Behaviours similar to those in figure 10 may be observed. However, as mentioned in § 4.2, in addition to $Fn_c^{(2i-1)}$, $Fn_c^{(2i)}$ will also affect the waves in these asymmetric cases about $y=0$. As a result, more sudden changes can be observed in figure 11.

Figure 10. Hydrodynamic forces on the cylinder when the ice sheet is under the symmetric edge conditions: (a) resistance; (b) lift. See figure 5 for legend.

Figure 11. Hydrodynamic forces on the cylinder when the ice sheet is under the asymmetric edge conditions: (a) resistance; (b) lift.

Computations are then carried out to investigate the effect of the channel width $b$. A comparison with the two-dimensional cases given in figure 4 of Li et al. (Reference Li, Wu and Shi2019) is also conducted. It should be noted that apart from that there are channel width $b$ and edge conditions here, all the physical parameters are chosen as the same as those in figure 4 of Li et al. (Reference Li, Wu and Shi2019). The results of clamped–clamped edges are given in figure 12. When $b=2H$, it can be seen that the resistance and lift are significantly different from those in the two-dimensional case. As $b$ increases, the results tend to those without the tank wall. When $b=8H$, very good agreement can be observed between these two cases. It ought to be pointed out at large $b$, for example $b=8H$, the critical Froude numbers $Fn_c^{(i)}$ ($i\geq 2$) still have some effects. However, the effects are obvious only when $Fn$ is very close to $Fn_c^{(i)}$, or the sharp variation of the force is almost like a vertical line, although this line is not included in the figure. In figure 13, the resistance and lift on the cylinder with the free–free edges are provided. It can be observed that even at small $b$, the results are already close to the two-dimensional ones. This suggests that when doing model tests in an ice-covered tank, the free–free edges may better resemble an unbounded ocean.

Figure 12. Hydrodynamic forces on the cylinder in an ice-covered channel with clamped–clamped edges at different channel widths: (a) resistance; (b) lift.

4.4. Deflection of the ice sheet induced by a submerged horizontal cylinder

We may also investigate the deflection of the ice sheet $\eta$ due to a submerged horizontal cylinder. The position of the centre line of cylinder is the same as that in § 4.3. In figure 14, the hydroelastic waves along the $x$-axis when the ice sheet is under clamped–clamped edges are plotted. As can be observed, many features of the waves in figure 14 are similar to those in figure 7, but there are also some differences. In figure 14(a), we can see that the wave profile for $Fn< Fn_c^{(1)}$ is symmetric about $x=x_0$, which can be shown through (3.53)–(3.55), while this is opposite to that in figure 7(a). As $Fn$ increases, the ice deflection $\eta$ at $x=0$ reaches a large negative value when $Fn$ is close to $Fn_c^{(1)}$, then it will rapidly become a positive value after $Fn$ passes the critical point. Similar sudden changes can be also observed in figures 14(c) and 14(d) when $Fn$ is near $Fn_c^{(3)}$ and $Fn_c^{(5)}$. However, the deflection above the cylinder in figure 14(c) before the critical point is a crest, while it is a trough in figures 14(a) and 14(d). As discussed previously, the results from the linear theory might become qualitatively inaccurate near one of the critical Froude numbers and the nonlinearity and other effects may need to be taken into account. However, this is only within a very small range. Also the critical Froude numbers predicted from the linear theory are still valid.

In order to show the influence of edge conditions on the wave profiles, typical wave patterns under six different edge constraints when $Fn_c^{(3)}< Fn< Fn_c^{(4)}$ are presented in figure 15. For waves under symmetric edges given in figures 15(a), 15(c) and 15(e), there is only one component in the longitudinal wave at the downstream region, while the waves at the upstream region contain two components. It can be found that the wave patterns under clamped–clamped and simply supported–simply supported edges are relatively similar, but the wave pattern under free–free is quite different. In figure 15(e), we can see that the downstream wave under free–free edges is not as obvious as that under the other two symmetric edges. Similar phenomena can be also observed from the transverse waves in figures 15(a), 15(c) and 15(e). From these results, we may infer that an ice-covered channel with free–free edges is closer to an unbounded region when compared with those with other edge constraints. For waves under asymmetric edge constraints shown in figures 15(b), 15(d) and 15(f). The longitudinal waves at the downstream region consist of two components, while the waves at the upstream region contains three components. Compared with the wave profiles under clamped–simply supported edges shown in figure 15(b), the asymmetry of the transverse waves is more obvious in cases with a free edge given in figures 15(d) and 15(f).

Figure 13. Hydrodynamic forces on the cylinder in an ice-covered channel with free–free edges at different channel widths: (a) resistance; (b) lift.

Figure 14. Ice sheet deflection $\eta$ along $x$-axis due to a submerged horizontal cylinder when the ice sheet is under clamped–clamped edges: (a) $Fn< Fn_c^{(1)}$; (b) $Fn_c^{(1)}< Fn< Fn_0$; (c) $Fn_c^{(2)}< Fn< Fn_c^{(3)}$; (d) $Fn_c^{(4)}< Fn< Fn_c^{(5)}$; with $Fn_c^{(1)}=0.8180$, $Fn_c^{(2)}=0.9981$, $Fn_c^{(3)}=1.3177$, $Fn_c^{(4)}=1.7682$, $Fn_c^{(5)}=2.3230$, $Fn_0=1.1469$.

Figure 15. Wave patterns under different edge conditions when $Fn_c^{(3)}< Fn< Fn_c^{(4)}$: (a) clamped–clamped edges, $Fn=1.5$ ($Fn_c^{(3)}=1.3177$ and $Fn_c^{(4)}=1.7682$); (b) clamped–simply supported edges, $Fn=1.4$ ($Fn_c^{(3)}=1.2569$ and $Fn_c^{(4)}=1.6824$); (c) simply supported–simply supported edges, $Fn=1.3$ ($Fn_c^{(3)}=1.2012$ and $Fn_c^{(4)}=1.6005$); (d) clamped-free edges, $Fn=1.3$ ($Fn_c^{(3)}=1.0962$ and $Fn_c^{(4)}=1.4454$); (e) free–free edges, $Fn=1.05$ ($Fn_c^{(3)}=0.9486$ and $Fn_c^{(4)}=1.1997$); (f) simply supported–free edges, $Fn=1.1$ ($Fn_c^{(3)}=1.0561$ and $Fn_c^{(4)}=1.3782$).