2.1 General approach based on the hodograph method

Usually, finding the function $w=w(z)$ directly is a complicated problem. Instead, Michell (Reference Michell1890) and Joukovskii (Reference Joukovskii1890) proposed the introduction of an auxiliary parameter plane, or $\unicode[STIX]{x1D701}$-plane, which was typically chosen as the upper half-plane. Then, they considered two functions, which were the complex potential $w$ and the function $\unicode[STIX]{x1D714}=-\text{ln}(\text{d}w/\text{d}z)$ in the parameter plane. When $w(\unicode[STIX]{x1D701})$ and $\unicode[STIX]{x1D714}(\unicode[STIX]{x1D701})$ are obtained as functions of $\unicode[STIX]{x1D701}$, the velocity field and the function mapping the parameter plane onto the physical plane can be determined as follows:

(2.4a,b)$$\begin{eqnarray}\frac{\text{d}w}{\text{d}z}=\exp [-\unicode[STIX]{x1D714}(\unicode[STIX]{x1D701})],\quad z(\unicode[STIX]{x1D701})=z_{0}+\int _{0}^{\unicode[STIX]{x1D701}}\left.\frac{\text{d}w}{\text{d}z}\right/\!\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701},\end{eqnarray}$$

where $z_{0}=z(0)$ in the physical plane corresponds to the origin in the $\unicode[STIX]{x1D701}$-plane.

The flow region beneath the free surface 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 double connected region becomes simply connected. As shown in figure 1(a), $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 and the free surface. The simply connected flow region including $O_{-}D_{+}$ and $O_{+}D_{-}$ is then transformed into the domain bounded by a rectangle $O_{-}O_{+}D_{-}D_{+}$ in the parameter plane. We choose the rectangle as the parameter domain to provide conformal mapping at points $O_{-}$, $O_{+}$, $D_{-}$, $D_{+}$. Indeed, the rectangle has the right angle at these points, the same as that at the corresponding points in the physical plane. For other forms of the parameter region, for which angles at points $O_{-}$, $O_{+}$, $D_{-}$, $D_{+}$ are different from a right angle, a singularity will appear, which makes it more difficult to find the solution. We may choose the coordinates of the rectangle vertices $O_{-},O_{+},D_{-},D_{+}$ as $(0,0)$, $(\unicode[STIX]{x03C0},0)$, $(\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$ and $(0,\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$, respectively, in the parameter plane shown in figure 1(b). Here, $\unicode[STIX]{x1D70F}$ is an imaginary number. The horizontal length of the rectangle is then equal to $\unicode[STIX]{x03C0}$ and its vertical length is equal to $\unicode[STIX]{x03C0}|\unicode[STIX]{x1D70F}|/4$.

In the flow region, there are two stagnation points marked $A$, where two streamlines merge into one, and $B$, where a streamline splits into two branches. Positions of these points in parameter plane $\unicode[STIX]{x1D701}$, $\unicode[STIX]{x1D701}=a$, $\unicode[STIX]{x1D701}=b$ as well as the position of point $C$, $\unicode[STIX]{x1D701}=c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$, corresponding to infinity in the physical plane, should be determined from additional conditions. The interval $0\leqslant \unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$ on the real axis corresponds to the body boundary. The interval $c<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$, $\unicode[STIX]{x1D702}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ corresponds to part of the free surface $D_{-}C_{-}$, and the interval $0\leqslant \unicode[STIX]{x1D709}<c$, $\unicode[STIX]{x1D702}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$, corresponds to the other part of the free surface $D_{+}C_{+}$.

It should be noticed that points $C_{-}$ as $x\rightarrow -\infty$ and $C_{+}$ as $x\rightarrow +\infty$ have been transformed to the same point $C$ in the parameter plane. As discussed in the first paragraph of this section, $y\rightarrow h$ at $C_{-}$ while it is a wavy function at $C_{+}$. Thus some treatments are needed to ensure they can become the same point in the parameter plane. For this reason, in the region between points $T_{1}$ and $T_{2}$ in figure 1(a), an artificial curved plate is placed on the free surface, which applies a pressure on the surface. The logic of this method is that, for steady free-surface flow past a body in infinite water depth, the body does not create a wave far upstream but it does far downstream. If the artificial body is put sufficiently behind the physical body, it does not affect the flow near the physical body significantly, as it does not generate a wave there. On the other hand, when its shape is adjusted properly, the wave generated by its presence can cancel the wave generated by the physical body far downstream, leading to a flat surface as $x\rightarrow \infty$. This is similar to what is used by Faltinsen & Semenov (Reference Faltinsen and Semenov2008). In Vanden-Broeck (Reference Vanden-Broeck2010) the fact that the free surface is infinite is reflected in an integral derived, where the integration limit tends to infinity. The limit is then truncated, which is equivalent to treating the free surface far downstream as a flat surface. To adopt this technique here, we can also choose the artificial body to be a flat plate. However, the plate has to be placed farther away than a curved body, which changes the free surface gradually rather than abruptly.

The flux of wave energy downstream is related to the force applied to the fluid (Wehausen & Liatone Reference Wehausen and Liatone1960; Newman Reference Newman1977). If the total force due to the submerged body and the artificial plate is zero, then there will be no wave at $x\rightarrow \infty$. As a result, $C_{+}$ and $C_{-}$ in the physical plane can be transformed to the same point $C$ in the parameter plane. The shape of the curved plate will be part of the solution. When it is applied, the velocity magnitude $v$ along $T_{1}T_{2}$ gradually tends to 1 towards $x\rightarrow x_{T2}$, and remains 1 beyond this, or $x_{T2}\leqslant x<\infty$. We should notice that the effect of the artificial plate on the flow near the real body will diminish as $x_{T1}$ increases. This is because the plate will not generate waves at its far upstream.

When solving free boundary problems, the shape of the auxiliary parameter region is usually chosen with the aim of obtaining the solution to the problem in the simplest form with the minimal number of singular points at which the transformation of the parameter region onto the complex potential $w$ and the function $\text{d}w/\text{d}z$ regions is not conformal. In the case of a simply connected domain, the fluid region is commonly converted to the first quadrant of the parameter plane or half-plane (e.g. Gurevich Reference Gurevich1965). In the case of the doubly connected flow region here, the additional corner points appear at the intersections of the two sides of the cut and the flow boundary. To correspond to these corner points, it is more convenient to choose the parameter plane in the form of a rectangle. On the other hand, when the integral equation is constructed in the half-plane or the first quadrant, polynomial functions are usually used. Here, for the rectangular domain, the polynomial functions will be replaced by Jacobi’s theta functions (Abramowitz & Stegun Reference Abramowitz and Stegun1964), which are quasi-doubly periodic functions. Jacobi’s functions have been used to solve free-surface problems involving doubly connected flow regions, for example by Birkhoff & Zarantonello (Reference Birkhoff and Zarantonello1957), Gurevich (Reference Gurevich1965) and Terentiev, Kirschner & Uhlman (Reference Terentiev, Kirschner and Uhlman2011).

2.2 Expression for complex velocity, $\text{d}w/\text{d}z$, defined in the $\unicode[STIX]{x1D701}$-plane

At this stage, we denote the angle of the velocity direction along the body as $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ and the velocity magnitude on the free surface as $v(\unicode[STIX]{x1D709})$. With these notations, we have the following boundary-value problem for the function of complex velocity, $\text{d}w/\text{d}z$:

(2.5)$$\begin{eqnarray}\displaystyle & \displaystyle \left|\frac{\text{d}w}{\text{d}z}\right|=v(\unicode[STIX]{x1D709}),\quad 0\leqslant \unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0},\,\unicode[STIX]{x1D702}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4, & \displaystyle\end{eqnarray}$$

(2.6)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D712}(\unicode[STIX]{x1D709})=\arg \left(\frac{\text{d}w}{\text{d}z}\right)=\left\{\begin{array}{@{}ll@{}}-\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709}), & 0\leqslant \unicode[STIX]{x1D709}<a,\,\unicode[STIX]{x1D702}=0,\\ -\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})-\unicode[STIX]{x03C0}, & a<\unicode[STIX]{x1D709}<b,\,\unicode[STIX]{x1D702}=0,\\ -\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})-2\unicode[STIX]{x03C0}, & b<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0},\,\unicode[STIX]{x1D702}=0,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.7)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}w}{\text{d}z}(\unicode[STIX]{x1D709}=0,\text{i}\unicode[STIX]{x1D702})=\frac{\text{d}w}{\text{d}z}(\unicode[STIX]{x1D709}=\unicode[STIX]{x03C0},\text{i}\unicode[STIX]{x1D702}),\quad 0\leqslant \text{i}\unicode[STIX]{x1D702}\leqslant \unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4. & \displaystyle\end{eqnarray}$$

In (2.6), the argument of the complex velocity has jumps equal to $-\unicode[STIX]{x03C0}$ at stagnation points $A$ ($\unicode[STIX]{x1D701}=a$) and $B$ ($\unicode[STIX]{x1D701}=b$) due to the jump in 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 of the rectangle.

The problem is to find the function $\text{d}w/\text{d}z$ in the parameter domain, which satisfies the given boundary conditions (2.5)–(2.7). Based on the Cauchy integral, it is possible to derive an expression for $\text{d}w/\text{d}z$ based on an integral along the boundary involving only the given condition. Here, we may follow a different approach. With the aim of applying the special points method (Gurevich Reference Gurevich1965), we decompose the body surface $O_{-}ABO_{+}$ in the parameter plane into $N$ intervals $E_{j-1},E_{j}$, with $j=1,\ldots ,N$, and assume that the argument of the velocity is piecewise constant, or $\unicode[STIX]{x1D712}(\unicode[STIX]{x1D709})\equiv \unicode[STIX]{x1D712}_{j-1}$ for $\unicode[STIX]{x1D709}\in (\unicode[STIX]{x1D709}_{j-1}\unicode[STIX]{x1D709}_{j})$, so that the velocity argument changes stepwise from $\unicode[STIX]{x1D712}_{j-1}$ to $\unicode[STIX]{x1D712}_{j}$ at vertex $E_{j}$. Here, $\unicode[STIX]{x1D709}_{0}=0$ and $\unicode[STIX]{x1D709}_{N}=\unicode[STIX]{x03C0}$ are defined. The order of the singularity of the function $\text{d}w/\text{d}z$ at point $E_{j}$ can be determined by analysing the behaviour of its argument when we are moving along an infinitesimal semi-circle centred at $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{j}$ in the counter-clockwise direction, as shown in figure 1(b). In such a way, the argument of ($\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j}$) changes by $\unicode[STIX]{x03C0}$, while the corresponding change in the argument of $\text{d}w/\text{d}z$ is $-(\unicode[STIX]{x1D712}_{j}-\unicode[STIX]{x1D712}_{j-1})=-\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}$. The theta function $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701})$ (see appendix A) to be used to derive $\text{d}w/\text{d}z$ has only one simple zero at point $\unicode[STIX]{x1D701}=0$ in the rectangle $O_{-}O_{+}D_{-}D_{+}$ (Whittaker & Watson Reference Whittaker and Watson1927), or $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701})\sim \unicode[STIX]{x1D701}$ near $\unicode[STIX]{x1D701}=0$. The argument of function $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})$ near point $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}_{j}$, $\arg [\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{j})]$, behaves in the same way as $\arg (\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})$. Therefore, near $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}_{j}$, if we write

(2.8)$$\begin{eqnarray}\frac{\text{d}w}{\text{d}z}\sim [\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})]^{-\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}/\unicode[STIX]{x03C0}},\quad j=1,\ldots ,N,\end{eqnarray}$$

then the change of the argument of $\text{d}w/\text{d}z$ is $-\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}$. Here (2.8) can be either zero or singular, depending on the sign of the exponent. To take into account the change of $\text{d}w/\text{d}z$ over each $\unicode[STIX]{x1D709}_{j}$, we may write

(2.9)$$\begin{eqnarray}\frac{\text{d}w}{\text{d}z}\sim \mathop{\prod }_{j=1}^{N}[\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})]^{-\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}/\unicode[STIX]{x03C0}},\quad \unicode[STIX]{x0394}\unicode[STIX]{x1D712}=\unicode[STIX]{x1D712}_{j}-\unicode[STIX]{x1D712}_{j-1},\,j=1,\ldots .N,\end{eqnarray}$$

when $z(\unicode[STIX]{x1D701})$ moves along the body surface. We also decompose the free surface $C_{-}D_{-}D_{+}C_{+}$ into $M$ intervals $F_{i-1}F_{i}$, with $i=1,\ldots ,M$, and assume that the velocity modulus within each interval $(\tilde{\unicode[STIX]{x1D709}}_{j-1},\tilde{\unicode[STIX]{x1D709}}_{j})$ is constant, $v(\tilde{\unicode[STIX]{x1D709}})\equiv v_{i}$. Similarly, the velocity modulus changes stepwise from $v_{i-1}$ to $v_{i}$ at point $F_{i}$. Here $\tilde{\unicode[STIX]{x1D709}}_{0}=\unicode[STIX]{x03C0}$ and $\tilde{\unicode[STIX]{x1D709}}_{M}=0$ are defined. On the basis of the above consideration, we consider the properties of the function $(\text{d}w/\text{d}z)^{i}=\exp [-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D701})]=\exp (\unicode[STIX]{x1D703}+\text{i}\ln v)$, where $\unicode[STIX]{x1D703}=-\text{arg}(\text{d}w/\text{d}z)$, and $\unicode[STIX]{x1D714}(\unicode[STIX]{x1D701})$ is defined in (2.4). The argument of this function is the logarithm of the velocity modulus, and its magnitude is the exponent of the velocity argument. Thus, moving along an infinitesimal semi-circle in the $\unicode[STIX]{x1D701}$-plane centred at $\unicode[STIX]{x1D701}_{i}=\tilde{\unicode[STIX]{x1D709}}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ in the counter-clockwise direction, $\arg [\unicode[STIX]{x1D701}-(\tilde{\unicode[STIX]{x1D709}}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)]$ changes by $\unicode[STIX]{x03C0}$, while the corresponding change in the argument of the function $(\text{d}w/\text{d}z)^{i}$ is $-(\ln v_{i}-\ln v_{i-1})$. Therefore, the function behaves near point $\unicode[STIX]{x1D701}_{i}=\unicode[STIX]{x1D709}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ as

(2.10)$$\begin{eqnarray}\left(\frac{\text{d}w}{\text{d}z}\right)^{i}\sim [\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{i}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)]^{-(1/\unicode[STIX]{x03C0})\ln (v_{i}/v_{i-1})},\quad i=1,\ldots ,M.\end{eqnarray}$$

Similarly to (2.9), we may then write

(2.11)$$\begin{eqnarray}\frac{\text{d}w}{\text{d}z}\sim \mathop{\prod }_{i=1}^{M}[\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{i}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)]^{(\text{i}/\unicode[STIX]{x03C0})\ln (v_{i}/v_{i-1})},\quad i=1,\ldots ,M.\end{eqnarray}$$

We now construct the expression for the complex velocity based on (2.9) and (2.11). Taking into account the requirement in (2.5) to (2.7), we may consider the following function

(2.12)$$\begin{eqnarray}f(\unicode[STIX]{x1D701})=\mathop{\prod }_{j=1}^{N}\left(\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}\right)^{-\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}/\unicode[STIX]{x03C0}}\mathop{\prod }_{i=1}^{M}\left(\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\tilde{\unicode[STIX]{x1D709}}_{i}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\tilde{\unicode[STIX]{x1D709}}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right)^{(\text{i}/\unicode[STIX]{x03C0})(v_{i}/v_{i-1})}.\end{eqnarray}$$

The denominator has no zeros or singularities in the rectangle $O_{-}O_{+}D_{-}D_{+}$ and its boundary, and therefore the zeros and singularities are those on the numerator, at $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}_{j}$, $j=1,\ldots ,N$ and $\unicode[STIX]{x1D701}=\tilde{\unicode[STIX]{x1D709}}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$, $i=1,\ldots ,M$, respectively.

As $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}+\unicode[STIX]{x03C0})=-\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701})$ (see appendix A), equation (2.7) is satisfied. We let $\unicode[STIX]{x1D701}=\tilde{\unicode[STIX]{x1D709}}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ on the free surface. Noticing that the nominator and denominator of the first product in (2.12) are complex conjugates as $\unicode[STIX]{x1D70F}$ is an imaginary number, we have

(2.13)$$\begin{eqnarray}\left|\frac{\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-\unicode[STIX]{x1D709}_{j}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-\unicode[STIX]{x1D709}_{j}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right|=1.\end{eqnarray}$$

Similarly, for point $\unicode[STIX]{x1D701}=\tilde{\unicode[STIX]{x1D709}}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ on the line $D_{-}D_{+}$, which is the mirror image of the point $\unicode[STIX]{x1D701}=\tilde{\unicode[STIX]{x1D709}}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ about $O_{-}O_{+}$ (see figure 1b), the second product in (2.12) can be evaluated as

(2.14)$$\begin{eqnarray}\displaystyle v_{0}\mathop{\prod }_{i=1}^{M}\left|\left(\frac{\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-\tilde{\unicode[STIX]{x1D709}}_{i})}{\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-\tilde{\unicode[STIX]{x1D709}}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}\right)^{(\text{i}/\unicode[STIX]{x03C0})\ln (v_{i}/v_{i-1})}\right| & = & \displaystyle v_{0}\exp \left\{\mathop{\sum }_{1}^{M}\ln \frac{v_{i}}{v_{i-1}}H(\tilde{\unicode[STIX]{x1D709}}_{i}-\tilde{\unicode[STIX]{x1D709}})\right\}\nonumber\\ \displaystyle & = & \displaystyle v_{0}\exp \left\{\mathop{\sum }_{1}^{k}\ln \frac{v_{i}}{v_{i-1}}\right\}=v_{k},\end{eqnarray}$$

where $H(\tilde{\unicode[STIX]{x1D709}}_{i}-\tilde{\unicode[STIX]{x1D709}})$ is the Heaviside step function. We can then verify that boundary conditions in (2.5) and (2.6) are satisfied by $f(\unicode[STIX]{x1D701})$ apart from a constant. Thus, the complex velocity may be written as

(2.15)$$\begin{eqnarray}\frac{\text{d}w}{\text{d}z}=C\mathop{\prod }_{j=1}^{N}\left(\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}\right)^{-\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}/\unicode[STIX]{x03C0}}\mathop{\prod }_{i=1}^{M}\left(\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\tilde{\unicode[STIX]{x1D709}}_{i}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\tilde{\unicode[STIX]{x1D709}}_{i}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right)^{(\text{i}/\unicode[STIX]{x03C0})(v_{i}/v_{i-1})},\end{eqnarray}$$

in which $C$ is a complex constant.

The products in (2.15) can be recast in the exponential form. We have

(2.16)$$\begin{eqnarray}\displaystyle \frac{\text{d}w}{\text{d}z} & = & \displaystyle C\exp \left[-\frac{1}{\unicode[STIX]{x03C0}}\mathop{\sum }_{1}^{N}\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j})}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}_{j}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\text{i}}{\unicode[STIX]{x03C0}}\mathop{\sum }_{1}^{M}\unicode[STIX]{x0394}\ln v_{i}\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\tilde{\unicode[STIX]{x1D709}}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\tilde{\unicode[STIX]{x1D709}}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}+\text{i}\unicode[STIX]{x1D712}(\unicode[STIX]{x03C0})\right].\end{eqnarray}$$

By letting $\unicode[STIX]{x0394}\unicode[STIX]{x1D712}_{j}=(\text{d}\unicode[STIX]{x1D712}/\text{d}\unicode[STIX]{x1D709})/\unicode[STIX]{x0394}\unicode[STIX]{x1D709}_{j}$ and $\unicode[STIX]{x0394}\ln v_{i}=\ln (v_{i}/v_{i-1})=(\text{d}\ln v/\text{d}\unicode[STIX]{x1D709})\unicode[STIX]{x0394}\unicode[STIX]{x1D709}_{i}$, and taking the limit as the step $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}$ in the argument and modulus of the velocity tends to zero and, correspondingly, $N\rightarrow \infty$ and $M\rightarrow \infty$, we obtain the integral formula

(2.17)$$\begin{eqnarray}\displaystyle \frac{\text{d}w}{\text{d}z} & = & \displaystyle v(\unicode[STIX]{x03C0})\exp \left[-\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D712}}{\text{d}\unicode[STIX]{x1D709}}\ln \left(\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}\right)\text{d}\unicode[STIX]{x1D709}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\text{i}}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}}\ln \left(\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right)\text{d}\unicode[STIX]{x1D709}+\text{i}\unicode[STIX]{x1D712}(\unicode[STIX]{x03C0})\right],\end{eqnarray}$$

which gives an expression for the mixed boundary-value problem (2.5)–(2.7) in the rectangle $O_{-}O_{+}D_{-}D_{+}$ based on the argument, $\arg [\text{d}w/\text{d}z]_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}=0}=\unicode[STIX]{x1D712}(\unicode[STIX]{x1D709})$, $0\leqslant \unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$, $\unicode[STIX]{x1D702}=0$ and the modulus, $|\text{d}w/\text{d}z|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}=v(\unicode[STIX]{x1D709})$. The important thing to notice here is that the right-hand side involves only those boundary values prescribed, either modulus or argument, and not the unknown part. Therefore, $\text{d}w/\text{d}z$ can be obtained directly once these prescribed values are known. The complex constant $C$ in (2.17) can be determined using the value of the argument of the function $\text{d}w/\text{d}z$ at point $O_{+}$ ($\unicode[STIX]{x1D701}=\unicode[STIX]{x03C0}$) and its modulus at point $D_{-}$ ($\unicode[STIX]{x1D701}=\unicode[STIX]{x03C0}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$). We note, that the integral formula (2.17) derived for the complex velocity, $\text{d}w/\text{d}z$, can be used to solve boundary-value problem (2.5)–(2.7) for an arbitrary complex function. By substituting the boundary conditions (2.5) and (2.6) into (2.17) and evaluating the first integral over the step change in the function $\unicode[STIX]{x1D712}(\unicode[STIX]{x1D709})$ at points $\unicode[STIX]{x1D709}=a$ and $\unicode[STIX]{x1D709}=b$, we obtain the expression for the complex velocity in the rectangle $O_{-}O_{+}D_{-}D_{+}$

(2.18)$$\begin{eqnarray}\displaystyle \frac{\text{d}w}{\text{d}z} & = & \displaystyle v_{D}\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-a)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-b)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-a)\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-b)}\exp \left[\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}}\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}\,\text{d}\unicode[STIX]{x1D709}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\text{i}}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}}\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\,\text{d}\unicode[STIX]{x1D709}-\text{i}\unicode[STIX]{x1D6FD}_{O}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{O_{-}}$ is the angle at point $O_{-}$ which is zero if point $O_{-}$ is the highest point of the body. Here, the relations between the theta functions $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709})$ and $\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709})$ have been used (see appendix A). The constant $v_{D}$ or the velocity magnitude at point $D_{+}$, is determined by satisfying the velocity at infinity, $\unicode[STIX]{x1D701}=c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$, which is 1, as it has been chosen as the reference velocity, or

(2.19)$$\begin{eqnarray}\left|\frac{\text{d}w}{\text{d}z}\right|_{\unicode[STIX]{x1D701}=c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}=1.\end{eqnarray}$$

2.3 Expression for derivative of the complex potential, $\text{d}w/\text{d}\unicode[STIX]{x1D701}$

For steady flows, the streamfunction $\unicode[STIX]{x1D713}=\text{Im}(w)$ takes constant values along the body and the free surface. According to Chaplygin’s special point method (see Gurevich (Reference Gurevich1965), chapter 1(5)), to determine the function $w=w(\unicode[STIX]{x1D701})$, it is sufficient to analyse all special points where the mapping is not conformal. These are stagnation points$A$ ($\unicode[STIX]{x1D701}=a$) and $B$ ($\unicode[STIX]{x1D701}=b$) and point $C$ ($\unicode[STIX]{x1D701}=c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$) corresponding to infinity in $w$-plane. The order of $w(\unicode[STIX]{x1D701})$ at these points can be determined by analysing the behaviour of the argument of $w(\unicode[STIX]{x1D701})$ in the vicinity of these points.

Moving in the counter-clockwise direction around an infinitesimal semi-circle centred at the point $\unicode[STIX]{x1D701}=a$, $\arg (\unicode[STIX]{x1D701}-a)$ changes by $\unicode[STIX]{x03C0}$ while the corresponding change of $\arg (w-w_{A})$ in the $w$-plane is $2\unicode[STIX]{x03C0}$. This means that the function $w(\unicode[STIX]{x1D701})$ behaves as

(2.20)$$\begin{eqnarray}w-w_{A}\sim \unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D701}-a)\end{eqnarray}$$

near point $A$. Similarly, at the stagnation point $B$ ($\unicode[STIX]{x1D701}=b$), the function $w(\unicode[STIX]{x1D701})$ behaves as

(2.21)$$\begin{eqnarray}w-w_{B}\sim \unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D701}-b).\end{eqnarray}$$

Moving in the counter-clockwise direction around an infinitesimal semi-circle centred at the point $\unicode[STIX]{x1D701}=c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$, $\arg (\unicode[STIX]{x1D701}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$ increases by $\unicode[STIX]{x03C0}$, while $\arg (w-w_{C})$ of the semi-circle of an infinite radius linking points $C_{+}$ and $C_{-}$ in the $w$-plane decreases by $\unicode[STIX]{x03C0}$. Therefore, the function $w(\unicode[STIX]{x1D701})$ behaves as

(2.22)$$\begin{eqnarray}w-w_{C}\sim \unicode[STIX]{x1D703}_{1}^{-1}(\unicode[STIX]{x1D701}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4).\end{eqnarray}$$

According to the symmetry principle, to satisfy the boundary condition $\text{Im}(w)=\text{const.}$ on the flow boundary, it is necessary to put special points of the same order at points $A^{\prime }$ ($\unicode[STIX]{x1D701}=a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2$) and $B^{\prime }$ ($\unicode[STIX]{x1D701}=b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2$), which are respectively symmetry points of $A$ and $B$ with respect to the side $D_{-}D_{+}$, and at point $C^{\prime }$ ($\unicode[STIX]{x1D701}=c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$) which is the symmetry point of $C$ with respect to side $O_{-}O_{+}$.

By differentiating the function $w(\unicode[STIX]{x1D701})$, we obtain the derivative $\text{d}w/\text{d}\unicode[STIX]{x1D701}$ which has the following zeros and poles:

(2.23)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D701}-a),\quad \unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D701}-b),\quad \unicode[STIX]{x1D703}_{1}^{-2}(\unicode[STIX]{x1D701}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4),\\ \unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D701}-a-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2),\quad \unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D701}-b-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2),\quad \unicode[STIX]{x1D703}_{1}^{-2}(\unicode[STIX]{x1D701}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\end{array}\right\}\end{eqnarray}$$

on boundaries of the rectangle and their symmetry lines. Similar to $\text{d}w/\text{d}z$, we may construct $\text{d}w/\text{d}\unicode[STIX]{x1D701}$ in the following form:

(2.24)$$\begin{eqnarray}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}=K^{\prime }\text{e}^{-\text{i}2(\unicode[STIX]{x1D701}-c)}\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-a)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-a-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-b)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-b-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}{\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D701}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D701}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)},\end{eqnarray}$$

where $K^{\prime }$ is a real scale factor. The term $\text{e}^{-\text{i}2(\unicode[STIX]{x1D701}-c)}$ is introduced into the above equation to satisfy kinematic boundary condition on the body and the free surface. In fact, on the free surface, $\unicode[STIX]{x1D701}=\tilde{\unicode[STIX]{x1D709}}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$, equation (2.24) takes the form

(2.25)$$\begin{eqnarray}\displaystyle \frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}} & = & \displaystyle K^{\prime }q^{-1/2}\text{e}^{-\text{i}2(\tilde{\unicode[STIX]{x1D709}}-c)}\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-a-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-b-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}^{2}(\tilde{\unicode[STIX]{x1D709}}-c)\unicode[STIX]{x1D717}_{1}^{2}(\tilde{\unicode[STIX]{x1D709}}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We notice that $\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-\tilde{\unicode[STIX]{x1D709}}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$ and $\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-\tilde{\unicode[STIX]{x1D709}}^{\prime }-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$ are complex conjugates to each other and

(2.26)$$\begin{eqnarray}\unicode[STIX]{x1D717}_{1}(\tilde{\unicode[STIX]{x1D709}}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)=\text{i}q^{-1/4}\text{e}^{-\text{i}(\tilde{\unicode[STIX]{x1D709}}-c)}\unicode[STIX]{x1D717}_{4}(\tilde{\unicode[STIX]{x1D709}}-c).\end{eqnarray}$$

Thus,

(2.27)$$\begin{eqnarray}\text{Im}(\text{d}w)=\text{Im}[(\text{d}w/\text{d}\unicode[STIX]{x1D701})_{\unicode[STIX]{x1D701}=\tilde{\unicode[STIX]{x1D709}}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}\,\text{d}\tilde{\unicode[STIX]{x1D709}}]=0.\end{eqnarray}$$

On the body surface, $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}$, we may use

(2.28)$$\begin{eqnarray}\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/2)=-\text{i}q^{1/4}\text{e}^{\text{i}\unicode[STIX]{x1D701}}\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701})\end{eqnarray}$$

to write (2.24) as

(2.29)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}=K\text{e}^{\text{i}(2c-a-b)}\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-a)\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709}-a)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-b)\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709}-b)}{\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D709}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D709}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}, & \displaystyle\end{eqnarray}$$

where $K=K^{\prime }q^{-1}$. The numerator is real following the definitions of the functions $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701})$ and $\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701})$. The denominator is also real since $\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D709}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$ and $\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D709}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)$ are complex conjugate quantities. Thus,

(2.30)$$\begin{eqnarray}\text{Im}(\text{d}w)=\text{Im}[(\text{d}w/\text{d}\unicode[STIX]{x1D701})_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}]=0,\end{eqnarray}$$

if $2c-a-b=0$ which will be shown later. Using this, equation (2.29) becomes

(2.31)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}=K\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-a)\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-a)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-b)\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-b)}{\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D701}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D701}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}. & \displaystyle\end{eqnarray}$$

Dividing (2.31) by (2.18), we obtain the derivative of the mapping function as

(2.32)$$\begin{eqnarray}\displaystyle \frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}} & = & \displaystyle \frac{K}{v_{D}}\frac{\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D701}-a)\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D701}-b)}{\unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D701}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D701}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\nonumber\\ \displaystyle & & \displaystyle \times \exp \left[-\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}}\ln \frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}\,\text{d}\unicode[STIX]{x1D709}\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{\text{i}}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}}\ln \frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\,\text{d}\unicode[STIX]{x1D709}+\text{i}\unicode[STIX]{x1D6FD}_{O}\right]\end{eqnarray}$$

whose integration along the intervals $0\leqslant \unicode[STIX]{x1D709}<c$ and $c<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$ at $\unicode[STIX]{x1D702}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$ in the $\unicode[STIX]{x1D701}$-plane provides the parts $D_{+}C_{+}$ and $D_{-}C_{-}$ of the free surface $C_{-}C_{+}$ in the $\unicode[STIX]{x1D701}$-plane, respectively. The parameters $a,b,c,\unicode[STIX]{x1D70F}$ and $K$, and the functions $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ and $\text{d}(\ln v)/\text{d}\unicode[STIX]{x1D709}$ 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.

2.4 System of equations for parameters $a,b,c,\unicode[STIX]{x1D70F}$ and $K$

At infinity, point $C_{-}C_{+}$ ($\unicode[STIX]{x1D701}_{C}=c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$), the velocity approaches unity (since this velocity is chosen as the reference velocity) and its direction is along the $x$-axis. Therefore, the argument of the complex velocity (2.18) at point $\unicode[STIX]{x1D701}_{C}$ should be equal to zero

(2.33)$$\begin{eqnarray}\arg \left(\frac{\text{d}w}{\text{d}z}\right)_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{C}}=1.\end{eqnarray}$$

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

(2.34)$$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}s_{b}}{\text{d}\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}=S_{b},\end{eqnarray}$$

where

(2.35)$$\begin{eqnarray}\frac{\text{d}s_{b}}{\text{d}\unicode[STIX]{x1D709}}=\left|\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}}\right|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}}.\end{eqnarray}$$

The free surface on the left- and right-hand sides at infinity has the same value of $y$-coordinate. This is also equivalent to the fact that the streamfunction $\unicode[STIX]{x1D713}=\text{Im}(w)$ is continuous across the cut, or $\text{Im}(w_{D_{-}})-\text{Im}(w_{D_{+}})=0$. By integrating $\text{Im}(\text{d}w/\text{d}\unicode[STIX]{x1D701})$ along $D_{-}D_{+}$ passing the point $\unicode[STIX]{x1D701}_{C}$ along a semi-circle $C^{\prime }$ of an infinitesimal radius $\unicode[STIX]{x1D700}$, at which $\text{d}w/\text{d}\unicode[STIX]{x1D701}$ in (2.31) has the second-order singularity, we have

(2.36)$$\begin{eqnarray}\displaystyle & & \displaystyle \text{Im}\left(\int _{\unicode[STIX]{x03C0}}^{c+\unicode[STIX]{x1D700}}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701}+\oint _{C^{\prime }}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701}+\int _{c-\unicode[STIX]{x1D700}}^{0}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\text{Im}\left(\oint _{C^{\prime }}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701}\right)=\text{Im}\left(\text{i}\unicode[STIX]{x03C0}\mathop{\text{Res}}_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{C}}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\text{Im}\left\{\text{i}\unicode[STIX]{x03C0}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D701}}\left[\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D701}_{C})^{2}\right]_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{C}}\right\}.\end{eqnarray}$$

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

(2.37)$$\begin{eqnarray}a+b=2c.\end{eqnarray}$$

The depth of submergence, $h$, and the flow rate, $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 flow rate $Q$, which is the integral of the derivative of the complex potential along the side $O_{-}D_{+}$

(2.38)$$\begin{eqnarray}\text{Im}\left(\int _{0}^{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701}\right)=Q.\end{eqnarray}$$

We may place a vortex with circulation $\unicode[STIX]{x1D6E4}$ at the centre of the cylinder, which can be non-dimensionalized as $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6E4}/(2\unicode[STIX]{x03C0}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, $\unicode[STIX]{x1D6FE}$ should be determined through the Kutta condition at the trailing edge.

Integrating $\text{d}w/\text{d}\unicode[STIX]{x1D701}$ along the body surface in the parameter plane, we have

(2.39)$$\begin{eqnarray}\text{Re}\left(\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701}\right)=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}.\end{eqnarray}$$

In the case $\unicode[STIX]{x1D6FE}\neq 0$, the real part of the potential, $\unicode[STIX]{x1D719}=\text{Re}(w)$, has a jump on the sides $O_{-}D_{-}$ and $O_{+}D_{+}$ of the cut, while the complex velocity, $\text{d}w/\text{d}z$ and the stream function $\unicode[STIX]{x1D713}=\text{Im}(w)$ are still continuous across the cut.

Equations (2.33)–(2.39) allow us to determine the unknown parameters $a,b,c,\unicode[STIX]{x1D70F}$ and $K$, which appear in the governing equations (2.18), (2.31) and (2.32), once the functions $v(\unicode[STIX]{x1D709})$ and $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ are specified.

2.5 Integro-differential equations for functions $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ and $v(\unicode[STIX]{x1D709})$

2.5.1 Kinematic boundary condition on the body surface

By integrating the modulus of the mapping function (2.32) 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

(2.40)$$\begin{eqnarray}s_{b}(\unicode[STIX]{x1D709})=\int _{0}^{\unicode[STIX]{x1D709}}\frac{\text{d}s_{b}}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\,\text{d}\unicode[STIX]{x1D709}^{\prime },\end{eqnarray}$$

where $\text{d}s_{b}/\text{d}\unicode[STIX]{x1D709}=|\text{d}z/\text{d}\unicode[STIX]{x1D701}|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}=0}$. Since the function $\unicode[STIX]{x1D6FD}_{b}(s_{b})$ is known, the function $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ can be determined from the following integro-differential equation:

(2.41)$$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}}=\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}s_{b}}\frac{\text{d}s_{b}}{\text{d}\unicode[STIX]{x1D709}}.\end{eqnarray}$$

By substituting $\text{d}z/\text{d}\unicode[STIX]{x1D701}$ from (2.32), this equation takes the form

(2.42)$$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}} & = & \displaystyle \unicode[STIX]{x1D705}[s_{b}(\unicode[STIX]{x1D709})]\frac{K}{v_{D}}\left|\frac{\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-a)\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-b)}{\unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D709}-c-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D709}-c+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right|\nonumber\\ \displaystyle & & \displaystyle \times \exp \left[-\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\ln \frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })}\,\text{d}\unicode[STIX]{x1D709}^{\prime }-\frac{\text{i}}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\ln \frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\,\text{d}\unicode[STIX]{x1D709}^{\prime }\right],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}(s_{b})=\text{d}\unicode[STIX]{x1D6FD}_{b}/\text{d}s_{b}$ is the curvature of the body.

2.5.2 Nonlinear boundary condition on the free surface: integral equation for the function $\text{d}\ln v/\text{d}\unicode[STIX]{x1D709}$

The velocity magnitude on the free surface can be obtained from the Bernoulli equation (2.3) in an iterative manner using the initial guess $y(\unicode[STIX]{x1D709})=h$ and

(2.43)$$\begin{eqnarray}v(\unicode[STIX]{x1D709})=\sqrt{1-\frac{2(y(\unicode[STIX]{x1D709})-h)}{F^{2}}}.\end{eqnarray}$$

Then, at each branch of the free surface, $C_{-}D_{-}$ ($c<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$) and $C_{+}D_{+}$ ($0\leqslant \unicode[STIX]{x1D709}<c$), its elevation can be obtained from

(2.44)$$\begin{eqnarray}y(\unicode[STIX]{x1D709})_{\{c<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0},0\leqslant \unicode[STIX]{x1D709}<c\}}=y_{D}+\text{Im}\left(\int _{\{\unicode[STIX]{x03C0},0\}}^{\unicode[STIX]{x1D709}}\left.\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}}\right|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}\,\text{d}\unicode[STIX]{x1D709}\right),\end{eqnarray}$$

where the integrand can be obtained from (2.32), and $y_{D}$ is the vertical coordinate of point $D$ and can be obtained from

(2.45)$$\begin{eqnarray}y_{D}=\text{Im}\left(\text{i}\int _{0}^{\unicode[STIX]{x03C0}|\unicode[STIX]{x1D70F}|/4}\left.\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}}\right|_{\unicode[STIX]{x1D701}=\text{i}\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D702}\right).\end{eqnarray}$$

This iteration converges well for the part $C_{-}D_{-}$ of the free surface for which $y(\unicode[STIX]{x1D709})$ changes monotonically, if $D$ is chosen appropriately. For the wavy part $C_{+}D_{+}$, convergence of the iteration is more difficult to achieve, which requires us to apply a more sophisticated approach.

Differentiating (2.3) with respect to the arc length coordinate along the free surface and taking into account that the slope of the free surface $\unicode[STIX]{x1D6FF}=\arcsin (\text{d}y/\text{d}s)=\unicode[STIX]{x03C0}+\unicode[STIX]{x1D6FD}$ is the angle between the unit tangential vector $\unicode[STIX]{x1D749}$ and the $x$-axis, we obtain

(2.46)$$\begin{eqnarray}F^{2}v^{2}\frac{\text{d}\ln v}{\text{d}s}-\sin \unicode[STIX]{x1D6FD}=0,\end{eqnarray}$$

where the velocity direction $\unicode[STIX]{x1D6FD}$ can be determined from (2.18) by putting $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$,

(2.47)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D709}) & = & \displaystyle \unicode[STIX]{x1D6FD}_{O}-\text{Im}\left\{\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right.\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\text{Im}\left[\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right]\text{d}\unicode[STIX]{x1D709}^{\prime }\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{1}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\ln \left|\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })}\right|\text{d}\unicode[STIX]{x1D709}^{\prime }\right\}.\end{eqnarray}$$

Using the relationships

(2.48)$$\begin{eqnarray}\frac{\text{d}v}{\text{d}s}=\left.v\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}}\right/\!\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D709}},\end{eqnarray}$$

and

(2.49)$$\begin{eqnarray}\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D709}}=\left|\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D701}}\right|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}=\frac{1}{v(\unicode[STIX]{x1D709})}\left|\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\right|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}=\frac{K}{v(\unicode[STIX]{x1D709})}\left|\frac{\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D709}-c)\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-c)}\right|.\end{eqnarray}$$

We can rewrite (2.46) with the following integral equation in terms of the function $\text{d}\ln v/\text{d}\unicode[STIX]{x1D709}$ over the interval $0\leqslant \unicode[STIX]{x1D709}<c$:

(2.50)$$\begin{eqnarray}\displaystyle & \displaystyle F^{2}v^{3}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}}+P_{1}(\unicode[STIX]{x1D709})\sin \left[\frac{1}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\ln \left|\frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime })}\right|\text{d}\unicode[STIX]{x1D709}^{\prime }+P_{2}(\unicode[STIX]{x1D709})\right]=0, & \displaystyle\end{eqnarray}$$

where

(2.51a,b)$$\begin{eqnarray}v(\unicode[STIX]{x1D709})=v_{D}\exp \left(\int _{0}^{\unicode[STIX]{x1D709}}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\,\text{d}\unicode[STIX]{x1D709}^{\prime }\right),\quad P_{1}(\unicode[STIX]{x1D709})=K\left|\frac{\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{1}^{2}(\unicode[STIX]{x1D709}-c)\unicode[STIX]{x1D703}_{4}^{2}(\unicode[STIX]{x1D709}-c)}\right|,\end{eqnarray}$$

(2.52)$$\begin{eqnarray}\displaystyle P_{2}(\unicode[STIX]{x1D709}) & = & \displaystyle \text{Im}\left\{\ln \frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}^{\prime }}\text{Im}\left\{\ln \frac{\unicode[STIX]{x1D703}_{1}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{4}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right\}\text{d}\unicode[STIX]{x1D709}^{\prime }-\unicode[STIX]{x1D6FD}_{O}.\end{eqnarray}$$

The velocity modulus function $v(\unicode[STIX]{x1D709})$ is determined by (2.43) for $c<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$ and by (2.50) for $0\leqslant \unicode[STIX]{x1D709}<c$.

2.6 Numerical approach

Discretization of the body and the free surfaces in the $\unicode[STIX]{x1D701}$-plane. In the discrete form, the solution is sought on a fixed set of points $\unicode[STIX]{x1D709}_{j}$, $j=1,\ldots ,N$ distributed along the side $O_{-}O_{+}$, $0\leqslant \unicode[STIX]{x1D709}_{j}\leqslant \unicode[STIX]{x03C0}$, $\unicode[STIX]{x1D702}=0$ and a fixed set of points $\tilde{\unicode[STIX]{x1D709}}_{i}$ distributed along the intervals $c+\unicode[STIX]{x1D700}\leqslant \tilde{\unicode[STIX]{x1D709}}_{i}\leqslant \unicode[STIX]{x03C0}$, $i=1,\ldots ,M_{1}$ and $0\leqslant \tilde{\unicode[STIX]{x1D709}}_{i}\leqslant c-\unicode[STIX]{x1D700}$, $i=M_{1}+1,\ldots ,M$ on the side $D_{-}D_{+}$, $\unicode[STIX]{x1D702}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4$. These intervals correspond to parts $D_{-}C_{-}$ and $D_{+}C_{+}$ of the free surface $C_{-}C_{+}$.

The value of $\unicode[STIX]{x1D700}$ is chosen to provide the length of the free surface $C_{+}D_{+}$ such that $x_{T1}>5\unicode[STIX]{x1D706}$, where $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}F^{2}$ is the wavelength predicted by a linear theory. The number of nodes on the body and the free surface are chosen in the range $N=100{-}300$ and $M=500{-}1000$, respectively, based on the requirement of convergence and accuracy of the solution. For all calculated examples, the difference between the results from the above ranges of $N$ and $M$ starts from the fourth figure only. The points $\unicode[STIX]{x1D709}_{j}$ are distributed in a way to provide a higher density of points $s_{j}=s_{b}(\unicode[STIX]{x1D709}_{j})$ near stagnation points $A$ ($\unicode[STIX]{x1D701}=a$) and $B$ ($\unicode[STIX]{x1D701}=b$). We have $\unicode[STIX]{x1D709}_{j}=a[1-\cos (\unicode[STIX]{x03C0}(j-1)/(N_{A}-1))]/2$, $j=1,\ldots ,N_{A}$, $\unicode[STIX]{x1D709}_{j}=a+(b-a)[1-\cos (\unicode[STIX]{x03C0}(j-N_{A})/(N_{B}-N_{A}))]/2$, $j=N_{A}+1,\ldots ,N_{B}$ and $\unicode[STIX]{x1D709}_{j}=b+(\unicode[STIX]{x03C0}-b)[1-\cos (\unicode[STIX]{x03C0}(j-N_{B})/(N-N_{B}))]/2$, $j=N_{B}+1,\ldots ,N$, where $N_{A}=N/4$ and $N_{B}=3N/4$.

The distribution of the points $\tilde{\unicode[STIX]{x1D709}}_{i}$ is chosen similarly to provide a higher density of points $s_{i}=s(\tilde{\unicode[STIX]{x1D709}}_{i})$ on the free boundary closer to the body.

Discrete form of the integro-differential equation for the function $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$. The integrals in the equations are evaluated using the linear interpolation of the functions $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ and $\ln v(\tilde{\unicode[STIX]{x1D709}})$ on the segments $(\unicode[STIX]{x1D709}_{j-1},\unicode[STIX]{x1D709}_{j})$ and $(\tilde{\unicode[STIX]{x1D709}}_{i-1},\tilde{\unicode[STIX]{x1D709}}_{i})$, respectively. Then, the integrals in (2.18) are evaluated numerically over each segment using trapezoidal rule.

The method of successive approximations is adopted for solving integro-differential equation (2.42) and integral equation (2.50). In the discrete form, equation (2.42) becomes

(2.53)$$\begin{eqnarray}\frac{(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}_{b})_{j}^{(k+1)}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D709}_{j}}=\frac{\unicode[STIX]{x1D6FD}_{b}[s_{b}^{(k)}(\unicode[STIX]{x1D709}_{j})]-\unicode[STIX]{x1D6FD}_{b}[s_{b}^{(k)}(\unicode[STIX]{x1D709}_{j-1})]}{\unicode[STIX]{x0394}\unicode[STIX]{x1D709}_{j}},\quad j=1,\ldots ,N,\end{eqnarray}$$

where the arc length along the body, $s_{b}^{(k)}(\unicode[STIX]{x1D709})$, is evaluated using (2.40) with $(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}_{bj}/\unicode[STIX]{x0394}\unicode[STIX]{x1D709}_{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}$.

Distribution of the velocity magnitude in the transition region $x_{T1}\leqslant x\leqslant x_{T2}$. In the $\unicode[STIX]{x1D701}$-plane, interval $(x_{T1},x_{T2})$ corresponds to the interval $(\tilde{\unicode[STIX]{x1D709}}_{T1},\tilde{\unicode[STIX]{x1D709}}_{T2})$. The integral equation (2.50) is solved over the interval $0\leqslant \tilde{\unicode[STIX]{x1D709}}<\tilde{\unicode[STIX]{x1D709}}_{T1}$ and, due to the radiation condition adopted in the present method, the velocity magnitude at point $T_{2}$ should be $v(x_{T2})=1$. Therefore, in the transition region, $x_{T1}\leqslant x\leqslant x_{T2}$ we use a one parameter family for the velocity distribution

(2.54)$$\begin{eqnarray}v(\overline{x})=\left\{\begin{array}{@{}ll@{}}v_{T1}+(1-v_{T1})\overline{x}+A_{v}\sin (2\unicode[STIX]{x03C0}\overline{x})\exp (-\overline{x}), & 0\leqslant \overline{x}\leqslant 1,\\ 1, & 1<\overline{x}<\infty ,\end{array}\right.\end{eqnarray}$$

where $\overline{x}=(x-x_{T1})/(x_{T2}-x_{T1})$ and $A_{v}$ is a parameter. The velocity distribution on the interval $(\tilde{\unicode[STIX]{x1D709}}_{T1},c)$ that is downstream of the point $T_{1}$ does not significantly affect the free surface near the body. This effect will diminish when $x_{T1}$ is sufficiently large. The parameter $A_{v}$ in (2.54) is chosen to further reduce this effect and obtain zero slope of the free surface at point $x_{M_{1}+1}$, or $\unicode[STIX]{x1D6FD}(\tilde{\unicode[STIX]{x1D709}}_{M_{1}+1})=0$.

The discrete form of the integral equation for the function $\text{d}\ln v/\text{d}\unicode[STIX]{x1D709}$. In the discrete form (2.46) can be written as follows:

(2.55)$$\begin{eqnarray}\displaystyle & \displaystyle F^{2}(\overline{v}_{i}^{3})^{(k)}\frac{\unicode[STIX]{x0394}(\ln v)_{i}^{(k+1)}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{i}}-\unicode[STIX]{x1D6FD}_{i}^{(k+1)}=\sin \unicode[STIX]{x1D6FD}_{i}^{(k)}-\unicode[STIX]{x1D6FD}_{i}^{(k)}. & \displaystyle\end{eqnarray}$$

Due to nonlinearity of (2.55), it is solved through an iteration procedure. The system of linear equations to determine the unknown $\unicode[STIX]{x0394}\ln v_{j}$ at the $(k+1)$th iteration can be written as follows:

(2.56)$$\begin{eqnarray}a_{ij}\unicode[STIX]{x0394}\ln v_{j}=b_{j},\quad i,j=1,\ldots ,n,\end{eqnarray}$$

where $n=M-M_{T}$ is the number of nodes $\tilde{\unicode[STIX]{x1D709}}_{i}$ between points $D_{+}$ and $T_{1}$, $\tilde{\unicode[STIX]{x1D709}}_{M_{T}}=\tilde{\unicode[STIX]{x1D709}}_{T1}$,

(2.57a-d)$$\begin{eqnarray}\displaystyle a_{ij}=A_{ij}+\unicode[STIX]{x1D6FF}_{ij}F^{2}(\overline{v}_{i}^{3})^{(k)}\frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{i}},\quad \overline{v}_{i}^{(k)}=\frac{v_{i-1}^{(k)}+v_{i}^{(k)}}{2},\quad \overline{\unicode[STIX]{x1D709}}_{i}=\frac{\tilde{\unicode[STIX]{x1D709}}_{i-1}+\tilde{\unicode[STIX]{x1D709}}_{i}}{2},\quad \unicode[STIX]{x1D6FF}_{ij}=\left\{\begin{array}{@{}ll@{}}1, & i=j,\\ 0, & i\neq j.\end{array}\right.\end{eqnarray}$$

(2.58)$$\begin{eqnarray}\displaystyle & \displaystyle A_{ij}=\frac{1}{\unicode[STIX]{x03C0}\unicode[STIX]{x0394}\tilde{\unicode[STIX]{x1D709}}_{j}}\int _{\tilde{\unicode[STIX]{x1D709}}_{j-1}}^{\tilde{\unicode[STIX]{x1D709}}_{j}}\ln \left|\frac{\unicode[STIX]{x1D717}_{1}(\overline{\unicode[STIX]{x1D709}}_{i}-\tilde{\unicode[STIX]{x1D709}}^{\prime })}{\unicode[STIX]{x1D717}_{4}(\overline{\unicode[STIX]{x1D709}}_{i}-\tilde{\unicode[STIX]{x1D709}}^{\prime })}\right|\text{d}\tilde{\unicode[STIX]{x1D709}}^{\prime },\quad \unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{i}=\text{Re}\left(\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\right)_{\unicode[STIX]{x1D701}=\overline{\unicode[STIX]{x1D709}}_{i}}\unicode[STIX]{x0394}\tilde{\unicode[STIX]{x1D709}}_{i}, & \displaystyle\end{eqnarray}$$

(2.59)$$\begin{eqnarray}\displaystyle & \displaystyle b_{i}=\sin \unicode[STIX]{x1D6FD}_{i}^{(k)}+\mathop{\sum }_{1}^{M_{T}}A_{ij}\unicode[STIX]{x0394}\ln v_{j}^{(k)},\quad \unicode[STIX]{x1D6FD}_{i}^{(k)}=-\mathop{\sum }_{1}^{M}A_{ij}\unicode[STIX]{x0394}\ln v_{j}^{(k)}+B_{i}, & \displaystyle\end{eqnarray}$$

(2.60)$$\begin{eqnarray}\displaystyle \displaystyle B_{i} & = & \displaystyle \unicode[STIX]{x1D6FD}_{O_{-}}-\text{Im}\left\{\ln \frac{\unicode[STIX]{x1D717}_{1}(\overline{\unicode[STIX]{x1D709}}_{i}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{1}(\overline{\unicode[STIX]{x1D709}}_{i}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(\overline{\unicode[STIX]{x1D709}}_{i}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{4}(\overline{\unicode[STIX]{x1D709}}_{i}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right\}\nonumber\\ \displaystyle & & \displaystyle -\,\mathop{\sum }_{1}^{N}\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}_{j}\text{Im}\left\{\frac{1}{\unicode[STIX]{x03C0}\unicode[STIX]{x0394}\tilde{\unicode[STIX]{x1D709}}_{j}}\int _{\tilde{\unicode[STIX]{x1D709}}_{j-1}}^{\tilde{\unicode[STIX]{x1D709}}_{j}}\ln \left|\frac{\unicode[STIX]{x1D717}_{1}(\overline{\unicode[STIX]{x1D709}}_{i}-\tilde{\unicode[STIX]{x1D709}}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(\overline{\unicode[STIX]{x1D709}}_{i}-\tilde{\unicode[STIX]{x1D709}}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right|\right\}.\end{eqnarray}$$

The iteration procedure is constructed as follows.

1. (i) Functions $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ and $v(\unicode[STIX]{x1D709})$ are initialized as $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})=2\unicode[STIX]{x1D709}$, $0\leqslant \unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0}$ and $v(\unicode[STIX]{x1D709})\equiv 1$, together with parameter $A_{v}=0$.

2. (ii) Equation (2.53) is solved by iteration for a given function $v(\unicode[STIX]{x1D709})$. For each iteration, equations (2.33)–(2.39) are solved to update $a$, $b$, $c$, $K$ and $\unicode[STIX]{x1D70F}$, respectively.

3. (iii) Apply external iteration to update $v(\unicode[STIX]{x1D709})$ by solving the system of (2.56).

4. (iv) Return to (ii) and continue until the free-surface shape downstream has converged.

5. (v) Update the velocity distribution upstream using (2.43) and (2.44) and return to (ii) until the velocity magnitude along the whole free surface has converged.

3.1 Flow past a submerged circular cylinder without gravity

The solution procedure is applied to solve the free-surface flow past a circular cylinder beneath the free surface without the effect of gravity. In such a case $F=\infty$ and (2.3) can be simplified as $v(\unicode[STIX]{x1D709})\equiv 1$. As the magnitude of the velocity is known on the free surface, only integro-differential equation (2.42), together with nonlinear equations (2.33)–(2.39) has to be solved to determine the function $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$. Figure 2 shows the streamline patterns at different circulations $\unicode[STIX]{x1D6FE}$ due to a vortex located at the centre of the cylinder. Due to the symmetry of the cylinder, these lines are symmetric about $x=0$.

Figure 2. Streamline patterns for a submerged circular cylinder at the flow rate $Q=1$ for different circulations: (a) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=0$; (b) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=3$; (c) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=-3$; (d) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=7.5$; (e) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=-7.5$.

The zero streamline is the line that passes through the stagnation point. It therefore splits into two lines at point $B$, which merge into one at point $A$. These two stagnation points are clearly shown in the figure. We may investigate the behaviour of the free surface at infinity using (2.44) in which $\text{Im}(\text{d}z/\text{d}\unicode[STIX]{x1D701})_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4}=\text{d}s/\text{d}\unicode[STIX]{x1D709}\sin \unicode[STIX]{x1D6FD}$

(3.1)$$\begin{eqnarray}y(\unicode[STIX]{x1D709})=y_{D}+\int _{0}^{\unicode[STIX]{x1D709}}\sin \unicode[STIX]{x1D6FD}\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709},\quad \unicode[STIX]{x1D709}\rightarrow c,\end{eqnarray}$$

where $\sin \unicode[STIX]{x1D6FD}=\text{d}y/\text{d}s$, and $s$ is the arc length along the free surface. The behaviour of $\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D709})$ on the free surface can be determined evaluating the leading order as $\unicode[STIX]{x1D709}\rightarrow c$. By taking into account (2.33), $\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D709})_{\unicode[STIX]{x1D709}=c}=0$, from (2.47) for $\unicode[STIX]{x1D709}\rightarrow c$ we can obtain from the Taylor expansion

(3.2)$$\begin{eqnarray}\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D709})\sim B_{1}(\unicode[STIX]{x1D709}-c),\quad \unicode[STIX]{x1D709}\rightarrow c,\end{eqnarray}$$

where

(3.3)$$\begin{eqnarray}\displaystyle B_{1} & = & \displaystyle -\text{Im}\left\{\frac{\unicode[STIX]{x1D717}_{1}^{\prime }(c-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(c-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}+\frac{\unicode[STIX]{x1D717}_{1}^{\prime }(c-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}(c-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{\unicode[STIX]{x1D717}_{4}^{\prime }(c-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(c-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}-\frac{\unicode[STIX]{x1D717}_{4}^{\prime }(c-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(c-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\right\}\nonumber\\ \displaystyle & & \displaystyle -\,\text{Im}\left\{\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}}{\text{d}\unicode[STIX]{x1D709}^{\prime }}[\unicode[STIX]{x1D717}_{1}^{\prime }(c-\unicode[STIX]{x1D709}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)-\unicode[STIX]{x1D717}_{4}^{\prime }(c-\unicode[STIX]{x1D709}^{\prime }+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)]\,\text{d}\unicode[STIX]{x1D709}^{\prime }\right\}\end{eqnarray}$$

and the dash on the theta function means the derivative. The leading order of the derivative $\text{d}s/\text{d}\unicode[STIX]{x1D709}$ is obtained from (2.49) taking into account that the function $\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D709}-c)$ along the free surface has a simple zero at $\unicode[STIX]{x1D709}=c$

(3.4)$$\begin{eqnarray}\displaystyle \frac{\text{d}s}{\text{d}\unicode[STIX]{x1D709}} & = & \displaystyle \frac{1}{v(\unicode[STIX]{x1D709})}\left|K\frac{\unicode[STIX]{x1D717}_{4}^{2}(\unicode[STIX]{x1D709}-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{4}^{2}(\unicode[STIX]{x1D709}-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D709}-c)\unicode[STIX]{x1D717}_{4}^{2}(\unicode[STIX]{x1D709}-c)}\right|\nonumber\\ \displaystyle & {\sim} & \displaystyle \frac{K_{1}}{v(c)\unicode[STIX]{x1D717}_{1}^{2}(\unicode[STIX]{x1D709}-c)}\sim \frac{K_{1}}{(\unicode[STIX]{x1D709}-c)^{2}},\quad \unicode[STIX]{x1D709}\rightarrow c,\end{eqnarray}$$

where

(3.5)$$\begin{eqnarray}K_{1}=\left|K\frac{\unicode[STIX]{x1D717}_{4}^{2}(c-a+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)\unicode[STIX]{x1D717}_{4}^{2}(c-b+\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D703}_{4}^{2}(0)}\right|.\end{eqnarray}$$

By using (3.2) and (3.4), from (3.1) can be obtained

(3.6)$$\begin{eqnarray}y(\unicode[STIX]{x1D709})\sim C_{1}\ln |\unicode[STIX]{x1D709}-c|,\quad \unicode[STIX]{x1D709}\rightarrow c,\end{eqnarray}$$

where $C_{1}=B_{1}/K_{1}$, and from (3.4)

(3.7)$$\begin{eqnarray}s(\unicode[STIX]{x1D709})\sim -K_{1}/(\unicode[STIX]{x1D709}-c),\quad \unicode[STIX]{x1D709}\rightarrow c.\end{eqnarray}$$

By eliminating the parameter $\unicode[STIX]{x1D709}$ from (3.6) and (3.7), we obtain

(3.8)$$\begin{eqnarray}y(s)\sim -C_{1}\ln |s|,\quad |s|\rightarrow \infty .\end{eqnarray}$$

If the constant $C_{1}=B_{1}/K_{1}\neq 0$, the $y$-coordinate of the free surface tends to infinity for $|s|\rightarrow \infty$. However, the value of coefficient $B_{1}$ in (3.3) very much depends on the combination of the parameters $a,b,c$ determined from the system of (2.33), (2.37) and (2.39).

Streamline patterns for a submerged circular cylinder at the flow rate $Q=1$ for different circulations are shown in figure 2. The free-surface elevation in figure 2(a) varies slowly in both directions away from the body, although it may not be obviously graphically visible. Figures 2(b) and 2(d) show streamlines for $\unicode[STIX]{x1D6FE}=3/(2\unicode[STIX]{x03C0})$ and $\unicode[STIX]{x1D6FE}=7.5/(2\unicode[STIX]{x03C0})$, respectively. Both stagnation points move towards the lower part of the cylinder, and the free-surface elevation at infinity tends to $-\infty$. For the circulation $\unicode[STIX]{x1D6FE}=-3/(2\unicode[STIX]{x03C0})$ and $\unicode[STIX]{x1D6FE}=-7.5/(2\unicode[STIX]{x03C0})$, the stagnation points move towards each other on the upper part of the cylinder, and the $y$-coordinate of the free surface tends to $+\infty$ at infinity. The flow patterns for larger flow rate $Q$ are shown in figure 3, which corresponds to large $y(x)_{x=0}$, or the distance between the free surface right above the body and its centre. As expected, the effects of the body and circulation on the local free-surface elevation become smaller than for the case shown in figure 2.

Figure 3. The streamline patterns for submerged circular cylinder at the flow rate $Q=10$ for different circulations: (a) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=0$; (b) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=3$; (c) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=-3$; (d) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=7.5$; (e) $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}=-7.5$.

The force acting on the body can be obtained from the Blasius–Chaplygin theorem

(3.9)$$\begin{eqnarray}\overline{R}_{f}=X_{f}-\text{i}Y_{f}=\frac{\text{i}\unicode[STIX]{x1D70C}U^{2}R}{2}\oint _{C_{b}}\left(\frac{\text{d}w}{\text{d}z}\right)^{2}\,\text{d}z=-\frac{\text{i}\unicode[STIX]{x1D70C}U^{2}R}{2}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}w}{\text{d}z}\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D701}}\,\text{d}\unicode[STIX]{x1D701},\end{eqnarray}$$

where $C_{b}$ is a closed contour enclosing the body and the integration follows the positive (anti-clockwise) direction. The sign ‘$-$’ in the last term appears due to the fact that integration route in the parameter plane from point $O_{-}$ to $O_{+}$ is in the opposite direction. Drag and lift coefficients are defined as follows:

(3.10a,b)$$\begin{eqnarray}C_{D}=\frac{2X_{f}}{\unicode[STIX]{x1D70C}U^{2}R},\quad C_{L}=\frac{2Y_{f}}{\unicode[STIX]{x1D70C}U^{2}R}.\end{eqnarray}$$

The dependences of the lift coefficient on circulation $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}$ at different $Q$ are shown in figure 4. For the large $Q$ corresponding to large $y(0)$, the lift coefficient $C_{L}/2\rightarrow -2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}$, which corresponds to the Kutta–Joukovskii theorem in the case without a free surface. For smaller $Q$, corresponding to smaller $y(0)$, the lift coefficient is smaller. Note, that the drag force is zero due to the symmetry of the flow about the $y$-axis.

Figure 4. Lift force versus circulation for different flow rates between the body and the free surface: $Q=0.5$ (dash-dotted line), $Q=1.0$ (dotted line), $Q=2.0$ (dashed line), $Q=20$ (solid line).

3.2 Flow past a submerged circular cylinder with gravity

For the case with gravity, we first reconsider those cases in Scullen & Tuck (Reference Scullen and Tuck1995), and their results obtained from the numerical method will be used for comparison.

The complete solution procedure presented in § 2.6 is used to predict the flow with gravity. For verification purposes we compare the obtained solution with the Scullen & Tuck (Reference Scullen and Tuck1995) results based on their numerical method. They represented the flow potential as a superposition of the potentials corresponding to the uniform flow, sources/sink and vortices placed on the surfaces external to the fluid domain. The singularities due to the sources were distributed at some distance over the free surface and within the circular cylinder close to its boundary. By using the dynamic and kinematic boundary conditions, the problem was reduced to a system of nonlinear equations solved by the Newton method. Scullen & Tuck (Reference Scullen and Tuck1995) presented results for various Froude numbers at depths of submergence $h\geqslant 4$.

Figure 5. Free-surface profiles for Froude number $F=2.75$ and depth of submergence $h=7.55$, or $F_{h}=1$, at different circulations: (a) $\unicode[STIX]{x1D6FE}=0$, (solid line), (b) $\unicode[STIX]{x1D6FE}=0.5\unicode[STIX]{x1D6FE}_{0}$ (dashed line), (c) $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{0}$ (dotted line). Symbols correspond to the Scullen & Tuck (Reference Scullen and Tuck1995) results and the lines to the present ones.

In figure 5 we show the free-surface profiles for the Froude number $F=2.75$ and depth of submergence $h=7.55$, which gives $F_{h}=1$. They are compared with those of Scullen & Tuck (Reference Scullen and Tuck1995), shown by the symbols. The circulation $\unicode[STIX]{x1D6FE}_{0}=\unicode[STIX]{x1D6E4}_{0}/(2\unicode[STIX]{x03C0}UR)=1/F^{2}$ is the one at which the solution from the first approximation theory is waveless far downstream or $x\rightarrow \infty$, and it predicts zero drag force (Tuck & Tulin Reference Tuck and Tulin1992). As can be seen from figure 5, the present results are in an excellent agreement with those of Scullen & Tuck (Reference Scullen and Tuck1995).

Figure 6. Free-surface profiles for different depths of submergence and Froude numbers: (a) $F=5$, (b) $F=3$, (c) $F=2$, (d) $F=1$.

Further calculations are made for different submergences and Froude numbers. The free-surface profiles are shown in figure 6, in which the $x$-coordinate is scaled by $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}F^{2}$, the wavelength from the linear theory. The centre of the cylinder is shown by the solid circle. For Froude numbers $F=5$ and $F=3$, shown in figures 6(a) and 6(b), respectively, the wave amplitude initially increases when $h$ decreases, reaches a maximum value at a submergence $h=h_{m}$ and then decreases as the cylinder further approaches the free surface. For $h<h_{m}$ the free-surface elevation above the cylinder is higher than the wave peak downstream. At very small depth, the downstream wave almost disappears. Note that, for $h<1$, the top of the cylinder is above the still water level. Due to inertia of the incoming liquid, the free surface rises near the cylinder, and the liquid between the zero streamline and the free surface passes the top of the cylinder. Then it moves down along the cylinder and merges with the main stream. Such flow may occur for those Froude numbers and submergences $h<1$, for which the free-surface elevation $y-h<F^{2}/2$, equivalent to the condition $v>0$ on the whole free surface.

For Froude number $F=2$ in figure 6(c), the free surface is less disturbed than that in the case of larger Froude number when $h>15$. However, at smaller submergence, or $h<6.2$, the wave height becomes larger. The curve for $h=3.93$ corresponds to the smallest depth of submergence for which the converged solution could be obtained from the iteration process. For the velocity magnitude equal to 0.5 at the wave crest, the wave elevation therefore is 75 % of the maximal possible elevation corresponding to $v=0$. As can be seen, for the case of $h=3.93$, the wave crest becomes sharper than the wave trough, and wave elevation departs from the sinusoidal curve, showing the important nonlinear effect. As the wave height increases, and the wavelength becomes smaller, which is consistent with nonlinear dispersion relationship (e.g. Lamb (Reference Lamb1932), art. 250). As the submergence further decreases, the wave steepness becomes larger. However, the results shown are for the cases within limit for which non-breaking waves can exist. Near or beyond the limit, the steady solution was not found, which was also observed and discussed by Scullen & Tuck (Reference Scullen and Tuck1995) and Faltinsen & Semenov (Reference Faltinsen and Semenov2008).

In our numerical method, when $h$ is nearer the critical value it becomes more difficult to obtain the solution, which is usually reflected by non-convergence during the iteration process. Convergence is obviously affected by the tolerance chosen for iteration. However, when the tolerance is sufficiently small, it no longer affects the limiting $h$ in the calculation. This was also observed and discussed by Haussling & Coleman (Reference Haussling and Coleman1979), Scullen & Tuck (Reference Scullen and Tuck1995) and by Faltinsen & Semenov (Reference Faltinsen and Semenov2008). In these cases, the wavelength and velocity continue to oscillate during the iterations, and the result does not converge.

For Froude $F=1$ in figure 6(d), the free surface is weakly disturbed for $h>6.2$. For $4.34<h<6.2$, the wave slope grows rapidly and for $h<4.34$ steady wave patterns were not obtained.

Figure 7. Wave height versus the depth of submergence for various Froude numbers: $F=5$ (solid line), $F=3$ (dashed line), $F=2$ (dotted line), $F=1$ (dot-dashed line).

The downstream wave height, $h_{w}=y_{cr}-y_{tr}$ versus submergence is shown in figure 7 for Froude numbers $F=5$, 3, 2 and 1. Here, $y_{cr}$ and $y_{tr}$ are the vertical coordinates of the wave crest and wave trough, respectively. For $F=5$ the maximal wave height is smaller than that for $F=3$. In the figure, at smaller Froude numbers $F=2$ and $F=1$, the wave height keeps increasing when the submergence decreases. Then, when the submergence further decreases, no converged solution could be obtained, as discussed above. This is principally owing to a larger value of $h_{w}/\unicode[STIX]{x1D706}_{n}$, where $\unicode[STIX]{x1D706}_{n}$ is the actual length of the nonlinear wave. In the cases of $F=5$ and 3, the wave height $h_{w}$ first increases and then reaches a peak when submergence decreases. Beyond the peak $h_{w}$ decreases with $h$. The solution can still be obtained even for $h<0$, for which the centre of the cylinder is already above the still water.

The pressure coefficient along the cylinder

(3.11)$$\begin{eqnarray}c_{p}=\frac{2(P-P_{\infty })}{\unicode[STIX]{x1D70C}U^{2}}=c_{pd}-\frac{2y}{F^{2}}\end{eqnarray}$$

in which $c_{pd}$ is the dynamic pressure and the second term is the hydrostatic pressure. The dynamic pressure coefficient as a function of the arc length coordinate $s=S/R$ is shown in figure 8 for Froude number $F=5$. The values $s/\unicode[STIX]{x03C0}=0$ and $s/\unicode[STIX]{x03C0}=2$ correspond to the same top point of the cylinder, and $s/\unicode[STIX]{x03C0}=1$ corresponds to the bottom point. For the case without free surface, the pressure coefficients at the top and bottom points are the same, and they are the troughs of the curve. As the cylinder approaches the free surface, or $h$ decreases, the trough corresponding to the bottom point goes up slightly, and the one corresponding to the top point moves up much more. Therefore, on excluding the buoyancy we can expect a downward vertical force whose magnitude increases as the cylinder approaches the free surface. For $h<1$, the top point of the cylinder is above the undisturbed free surface. The free surface upstream will rise near the cylinder and the pass over the top of the body, similar to that shown in figure 2(a). The peaks of the pressure coefficient $c_{pd}$ without the free surface in figure 8 correspond to the stagnation points, and the position is affected only slightly by the free surface.

Figure 8. Pressure distribution along the cylinder for $F=5$ at different submergences: no free surface (solid line), $h=2.87$ (dashed line), $h=1.48$ (dotted line), $h=0.86$ (dash-dotted line) and $h=0.47$ (dash-dot-dot line).

Figure 9. (a) Scaled drag and lift coefficients against $F_{h}=F/\sqrt{h}$ at different Froude numbers, $F$: the first approximation theory (solid line) and the nonlinear theory for $F=5$ (dashed line), $F=3$ (dotted line) and $F=2$ (dash-dotted line, ends at $F_{h}=1.2$); (b) scaled drag (upper curves) and lift (lower curves) coefficients as functions of the depth of submergence: nonlinear (solid lines and solid markers) and the first approximation theory (dashed lines and opened markers) for Froude $F=5$ (●, ○), $F=3$ (▴, ▵) and $F=2$ (▪, ▫).

Through the first approximation solution (Lamb Reference Lamb1932; Kochin, Kibel & Roze Reference Kochin, Kibel and Roze1964) the explicit equations for the resistance and lift force coefficients are obtained as

(3.12)$$\begin{eqnarray}C_{D}^{\prime }h^{2}=\frac{X_{f}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}gR^{2}}h^{2}=C_{D}\frac{h^{2}F^{2}}{2\unicode[STIX]{x03C0}}=4\unicode[STIX]{x03C0}\left(\frac{1}{F_{h}^{2}}-\unicode[STIX]{x1D6FE}h\right)\exp ^{-2/F_{h}^{2}},\end{eqnarray}$$

(3.13)$$\begin{eqnarray}\displaystyle C_{L}^{\prime }h^{2} & = & \displaystyle \frac{Y_{f}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}gR^{2}}h^{2}=C_{L}\frac{h^{2}F^{2}}{2\unicode[STIX]{x03C0}}=-\unicode[STIX]{x1D6FE}F_{h}^{2}h(2h^{2}+\unicode[STIX]{x1D6FE}h-1)+4\unicode[STIX]{x1D6FE}h-\frac{1}{2}F_{h}^{2}-1-\frac{2}{F_{h}^{2}}\nonumber\\ \displaystyle & & \displaystyle +\,4\left(\frac{1}{F_{h}^{2}}-\unicode[STIX]{x1D6FE}h\right)^{2}\exp ^{-2/F_{h}^{2}}Ei_{1}\left(\frac{2}{F_{h}^{2}}\right),\end{eqnarray}$$

where $Ei_{1}(z)$ is the exponential integral function. The coefficients $C_{D}^{\prime }$ and $C_{L}^{\prime }$ represent the force coefficients non-dimensionalized by the buoyancy of the cylinder, $\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}gR^{2}$. Results for the force at $\unicode[STIX]{x1D6FE}=0$ are obtained from the present fully nonlinear theory and are given in figure 9(a), in which $C^{\prime }h^{2}$ is plotted against submergence based on the Froude number, $F_{h}$, at different radius based $F$. For the first approximation, at $\unicode[STIX]{x1D6FE}=0$, $C^{\prime }h^{2}$ in (3.12) depends on $F_{h}=F/\sqrt{h}$ only. This means that, provided $F_{h}$ is the same, $C^{\prime }h^{2}$ does not change with $F$. Thus results from the first approximation theory at $F=5$, 3 and 2, equation (3.12) are shown by the same solid line in figure 9(a). However, based on the present nonlinear solution, the results are very much affected by $F$ even at the same $F_{h}$. The nonlinear solution can be found for all three $F$ when $F_{h}<1.2$, and the results are in a fairly good agreement with the first approximation solution. For $F_{h}>1.2$, no solution is obtained for the case of $F=2$. The results from $F=3$ and $F=5$, which correspond to a body of smaller radius relative to $h$, show some discrepancy with the first approximation solution. The difference is larger for the former, as the body is larger relative to the depth and therefore has a larger nonlinear effect.

The peaks of $C_{D}^{\prime }h^{2}$ in figure 9(a) occur around the same $F_{h}$ at different $F$. The curves for $C_{D}^{\prime }h^{2}$ and $C_{L}^{\prime }h^{2}$ versus the depth of submergence are shown in figure 9(b). The position is shifted to the left for large $F$ according to the relation $F_{h}=F/\sqrt{h}$. The normalized lift force coefficient, $C_{L}^{\prime }h^{2}$, becomes noticeable when $h$ becomes relatively small. The variation is faster for $F=5$ and then for $F=3$. For $F=2$ the curve ends at relatively larger submergence.

3.3 Flow past a submerged hydrofoil

The solution procedure derived in § 2 is also applicable for a body with a sharp corner point. In such a case, the Kutta–Joukovskii condition is to be imposed at the corner to ensure the velocity and the pressure are finite there. As a result, the circulation $\unicode[STIX]{x1D6FE}$ should be obtained from such a condition, rather than being chosen randomly. Its chord length is chosen as the characteristic length $L$ in the following computation.

We consider the hydrofoil NACA0012 which has a thickness of the trailing edge 2.5 % of the hydrofoil thickness. In the present computations, we neglected the thickness of the trailing edge by taking it to be zero. The inner angle at the trailing edge is assumed as $\unicode[STIX]{x1D707}$. In order to satisfy the Kutta condition, the stagnation point $A$ and the trailing edge of the foil should coincide. With this assumption, the slope of the hydrofoil can be written as follows:

(3.14)$$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D6FD}_{b}^{\ast }(\unicode[STIX]{x1D709}), & 0\leqslant \unicode[STIX]{x1D709}<a,\,\unicode[STIX]{x1D702}=0,\\ \unicode[STIX]{x1D6FD}_{b}^{\ast }(\unicode[STIX]{x1D709})-\unicode[STIX]{x03C0}+\unicode[STIX]{x1D707}, & a<\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x03C0},\,\unicode[STIX]{x1D702}=0,\end{array}\right.\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{b}^{\ast }(\unicode[STIX]{x1D709})$ is a continuous function changed from 0 to $2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D707}$. By substituting (3.14) into (2.18) and evaluating the first integral over the step change in the function $\unicode[STIX]{x1D6FD}_{b}(\unicode[STIX]{x1D709})$ at point $\unicode[STIX]{x1D709}=a$, we obtain the expression for the complex velocity for the hydrofoil

(3.15)$$\begin{eqnarray}\displaystyle \frac{\text{d}w}{\text{d}z} & = & \displaystyle v_{D}\left[\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-a)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-a)}\right]^{\unicode[STIX]{x1D707}/\unicode[STIX]{x03C0}}\frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-b)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-b)}\exp \left[\frac{1}{\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\frac{\text{d}\unicode[STIX]{x1D6FD}_{b}^{\ast }}{\text{d}\unicode[STIX]{x1D709}}\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709})}\,\text{d}\unicode[STIX]{x1D709}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{\text{i}}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x03C0}}^{0}\frac{\text{d}\ln v}{\text{d}\unicode[STIX]{x1D709}}\ln \frac{\unicode[STIX]{x1D717}_{1}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}{\unicode[STIX]{x1D717}_{4}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D709}-\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}/4)}\,\text{d}\unicode[STIX]{x1D709}-\text{i}\unicode[STIX]{x1D6FD}_{O}\right].\end{eqnarray}$$

It is seen from the above equation that the complex velocity at the trailing edge has a zero of order $\unicode[STIX]{x1D707}/\unicode[STIX]{x03C0}$. For the case $\unicode[STIX]{x1D707}=0$ which is assumed below, the first product in (3.15) disappears and the velocity at the trailing edge becomes finite. Point $A$ will still be chosen at the trailing edge. However, this may not be a stagnation point anymore. The circulation $\unicode[STIX]{x1D6FE}$ is obtained from (2.39). The parameters $a,b,c,K$ and $\unicode[STIX]{x1D70F}$ are obtained from the system of (2.33)–(2.38) and the equation

(3.16)$$\begin{eqnarray}s_{OA}=\int _{0}^{a}\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709},\end{eqnarray}$$

where $\text{d}s/\text{d}\unicode[STIX]{x1D709}=|\text{d}z/\text{d}\unicode[STIX]{x1D701}|_{\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}}$, and $\unicode[STIX]{x1D709}=0$ corresponds to the point where the slope of the upper part of the foil is equal to zero. The system of integral equations from (2.42) to (2.50) has to be slightly modified by replacing the expression for the complex velocity (2.18) with (3.15).

Figure 10. Free-surface profiles past a NACA0012 hydrofoil at angle of attack $\unicode[STIX]{x1D6FC}=5^{\circ }$, Froude number $F=0.567$ and depth of submergence (a) $h/L=1.256$ and (b) $h/L=1.0345$ measured from the mid-chord point. Solid lines, present solution; dashed lines, Faltinsen & Semenov (Reference Faltinsen and Semenov2008) analytical/numerical approach; dotted lines, Landrini et al. (Reference Landrini, Lugni and Bertram1999) nonlinear numerical approach; symbol ○, measurements from Duncan (Reference Duncan1983).

The free-surface elevations for the flow past the NACA0012 hydrofoil at the angle of attack $\unicode[STIX]{x1D6FC}=5^{\circ }$ and Froude number based on the chord length of the foil $F=0.567$ are shown in figure 10 and compared with the experimental data of Duncan (Reference Duncan1983) and with the results obtained through the fully nonlinear numerical approach by Landrini et al. (Reference Landrini, Lugni and Bertram1999) and the nonlinear solution by Faltinsen & Semenov (Reference Faltinsen and Semenov2008). Agreement between the present nonlinear solution and that by Faltinsen & Semenov (Reference Faltinsen and Semenov2008), as well as with that from the numerical simulation of Landrini et al. (Reference Landrini, Lugni and Bertram1999), is quite good. The agreement with the experimental data at larger submergence is also quite good. However, there is some discrepancy between the experimental data and all of the calculated data for the smaller depth of submergence (figure 10b).

Figure 11. Effect of submergence on the pressure distribution along the lower and upper sides of the hydrofoil NACA0012 for Froude number $F=2$ and angle of attack $\unicode[STIX]{x1D6FC}=5^{\circ }$: no free surface (solid line), $h=0.74$ (dashed line), $h=0.33$ (dotted line) and $h=0.10$ (dash-dotted line).

The dynamic pressure coefficient on the hydrofoil is plotted against $x$ for the Froude number $F=2$ and angle of attack $\unicode[STIX]{x1D6FC}=5^{\circ }$ in figure 11. The origin is located in the middle point of the foil between the leading and trailing edges. The pressure coefficient with peak $c_{pd}=1$ at the stagnation point near the leading edge, $x=x_{B}$ gradually decreases along the lower side of the hydrofoil. On the upper side, the pressure coefficient rapidly drops to its minimal value and then gradually increases with $x$. The effect of submergence is most pronounced on the upper side. It increases as the depth decreases, and even becomes larger than that on the lower side at the same $x$ on some parts of the hydrofoil.

Figure 12. Effect of submergence on the lift coefficient for the hydrofoil NACA0012 at the angle of attack $\unicode[STIX]{x1D6FC}=5^{\circ }$ for Froude numbers $F=10$, 5, 3, 2 and 1 (lines marked by square, circle, triangle, reversed triangle and diamond, respectively). The solid lines and symbols correspond to integration of the pressure around the hydrofoil, and the dashed lines and open symbols correspond to the Kutta–Joukowskii theorem relating the lift force and circulation around the hydrofoil.

The effect of the submergence on the lift coefficient is shown in figure 12. The solid lines and symbols correspond to the lift coefficients obtained through the integration of the pressure around the hydrofoil, while the dashed lines and unfilled symbols are based on the Kutta–Joukowskii theorem ($Y_{f}=-\unicode[STIX]{x1D70C}U\unicode[STIX]{x1D6E4}$) for the unbounded flow domain, which in the present notation is $C_{L}=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}$, where $\unicode[STIX]{x1D6FE}$ is obtained from (2.39). From the figure, it is seen that the lift decreases as the hydrofoil approaches the free surface, which is consistent with the increase of the pressure on the upper side of the hydrofoil in figure 11. At small depths of submergence, the smaller the Froude number is, the smaller the lift coefficient becomes, which may also become negative.

The wave resistance of a submerged body can also be determined by the wave energy flux generated downstream (Wehausen & Liatone Reference Wehausen and Liatone1960). For linear sinusoidal waves with height $h_{w}$, this may be simplified as (Newman Reference Newman1977)

(3.17)$$\begin{eqnarray}C_{D}=\frac{h_{w}^{2}}{16F^{2}},\end{eqnarray}$$

which does not depend on the body shape. In figure 13, the coefficients versus depth $h$ obtained using (3.9) and (3.17) are shown by the solid and dashed lines, respectively. As the depth decreases, the wave height becomes larger and $C_{D}$ reaches its maximal value $C_{D\,max}$ at submergence $h_{m}$. For a further decrease of the depth, $h<h_{m}$, the wave amplitude and the resistance, correspondingly, become smaller. For smaller Froude number the coefficient $C_{D\,max}$ becomes larger. However, at very small Froude number ($F<1$ in the presented calculations), the maximal wave resistance cannot be reached due to wave breaking, as can be seen in figure 13 for $F=0.567$. Thus, similar to a smooth body, the steady solution is not always possible for a hydrofoil either.

Figure 13. Effect of submergence on the wave resistance of the hydrofoil NACA0012 at the angle of attack $\unicode[STIX]{x1D6FC}=5^{\circ }$ for Froude numbers $F=10$, 5, 3, 2, 1 and 0.567. The solid and dashed lines correspond to (3.10) and (3.17), respectively.