3.1 Solution of $h(x,z)$ for given ${\mathcal{X}}(x,z)$

For given ${\mathcal{X}}$ , $h$ can be directly determined from the boundary-value problem consisting of (2.7), (2.8), (2.10) and the radiation condition by the use of a Fourier transform. As an illustration and for simplicity, we first consider the 2D problem (in the $x$ – $z$ plane) which was solved by Dagan & Miloh (Reference Dagan and Miloh1982) in the case of a point source. Their key solution procedure is included here to assist in understanding the analysis in the case of a body (cf. § 3.3).

We write the solution of $h(x,z)$ in the form

(3.1) $$\begin{eqnarray}\displaystyle h(x,z)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{\infty }\tilde{h}(k)\text{e}^{-\text{i}kx}\text{e}^{|k|z}\,\text{d}k,\quad z\leqslant 0, & & \displaystyle\end{eqnarray}$$

where $\tilde{h}(k)$ is the Fourier transform of $h(x,z=0)$ ,

(3.2) $$\begin{eqnarray}\displaystyle \tilde{h}(k)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{\infty }h(x,0)\text{e}^{\text{i}kx}\,\text{d}x. & & \displaystyle\end{eqnarray}$$

For clarity in later reference, we denote the Fourier transform of $h^{\ast }(x,z)$ by

(3.3) $$\begin{eqnarray}\displaystyle \tilde{h}^{\star }(k)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{\infty }h^{\ast }(x,0)\text{e}^{\text{i}kx}\,\text{d}x. & & \displaystyle\end{eqnarray}$$

Here, $h(x,z)$ in (3.1) satisfies Laplace’s equation (2.7) and the deep water condition (2.10). The unknown function $\tilde{h}(k)$ can be found from the free-surface boundary condition (2.8) together with the radiation condition.

Upon taking the Fourier transform of (2.8), we have

(3.4) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{8}\int _{-\infty }^{\infty }\unicode[STIX]{x1D735}h^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}h\boldsymbol{\cdot }\unicode[STIX]{x1D735}h)\text{e}^{\text{i}kx}\,\text{d}x+\frac{1}{8}\int _{-\infty }^{\infty }\unicode[STIX]{x1D735}h\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}h^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}h)\text{e}^{\text{i}kx}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{8}\int _{-\infty }^{\infty }\unicode[STIX]{x1D735}h\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}h^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}h)\text{e}^{\text{i}kx}\,\text{d}x+A(\unicode[STIX]{x1D70F},k)\tilde{h}(k)=\tilde{P}(k)\quad \text{on }z=0,\end{eqnarray}$$

where

(3.5) $$\begin{eqnarray}\displaystyle \tilde{P}(k)=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{\infty }P(x)\text{e}^{\text{i}kx}\,\text{d}x & & \displaystyle\end{eqnarray}$$

and

(3.6) $$\begin{eqnarray}\displaystyle A(\unicode[STIX]{x1D70F},k)\equiv -\unicode[STIX]{x1D70F}^{2}-2\unicode[STIX]{x1D70F}k+|k|-k^{2}. & & \displaystyle\end{eqnarray}$$

By making use of the Faltung theorem (e.g. Dagan & Miloh Reference Dagan, Miloh and Inui1980, Reference Dagan and Miloh1982; Bracewell Reference Bracewell1986) for the cubic terms on the left-hand side in (3.4), we obtain

(3.7) $$\begin{eqnarray}\displaystyle A(\unicode[STIX]{x1D70F},k)\tilde{h}(k)-\frac{1}{4}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }a(k,\unicode[STIX]{x1D708},\unicode[STIX]{x1D706})\tilde{h}(\unicode[STIX]{x1D708})\tilde{h}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})\tilde{h}^{\star }(k-\unicode[STIX]{x1D706})\,\text{d}\unicode[STIX]{x1D708}\,\text{d}\unicode[STIX]{x1D706}=\tilde{P}(k), & & \displaystyle\end{eqnarray}$$

where

(3.8) $$\begin{eqnarray}\displaystyle a(k,\unicode[STIX]{x1D708},\unicode[STIX]{x1D706}) & = & \displaystyle \frac{1}{\unicode[STIX]{x03C0}}\left\{(\unicode[STIX]{x1D70F}+\unicode[STIX]{x1D708})[|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}|+|k-\unicode[STIX]{x1D706}|](k-\unicode[STIX]{x1D706})[(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})(k-\unicode[STIX]{x1D706})\right.\nonumber\\ \displaystyle & & \displaystyle -\,|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}||k-\unicode[STIX]{x1D706}|]-(k-\unicode[STIX]{x1D706}-\unicode[STIX]{x1D70F})[2\unicode[STIX]{x1D70F}\unicode[STIX]{x1D706}^{2}|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}|-2\unicode[STIX]{x1D70F}\unicode[STIX]{x1D708}|\unicode[STIX]{x1D708}|(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})\nonumber\\ \displaystyle & & \displaystyle +\,|\unicode[STIX]{x1D708}|(\unicode[STIX]{x1D708}-|\unicode[STIX]{x1D708}|)(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})(|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}|+|\unicode[STIX]{x1D708}|)]-\left[(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})\unicode[STIX]{x1D708}^{2}(k-\unicode[STIX]{x1D706})\right.\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})\unicode[STIX]{x1D708}(k-\unicode[STIX]{x1D706})^{2}-|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}||\unicode[STIX]{x1D708}|(k-\unicode[STIX]{x1D706})\unicode[STIX]{x1D708}-(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})|\unicode[STIX]{x1D708}||k-\unicode[STIX]{x1D706}|\unicode[STIX]{x1D708}\nonumber\\ \displaystyle & & \displaystyle \left.\left.-\,(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708})|\unicode[STIX]{x1D708}|(k-\unicode[STIX]{x1D706})|k-\unicode[STIX]{x1D706}|+|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}|\unicode[STIX]{x1D708}^{2}|k-\unicode[STIX]{x1D706}|+\frac{1}{2}|\unicode[STIX]{x1D706}-\unicode[STIX]{x1D708}||\unicode[STIX]{x1D708}|(k-\unicode[STIX]{x1D706})^{2}\right]\right\}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

From (3.7), we obtain the solution of $\tilde{h}(k)$ , which can be expressed in a form similar to the linear solution,

(3.9) $$\begin{eqnarray}\displaystyle \tilde{h}(k)=\frac{\tilde{P}(k)}{A_{g}(\unicode[STIX]{x1D70F},k)}, & & \displaystyle\end{eqnarray}$$

with

(3.10) $$\begin{eqnarray}\displaystyle A_{g}(\unicode[STIX]{x1D70F},k)=A(\unicode[STIX]{x1D70F},k)+d(\unicode[STIX]{x1D70F},k,f), & & \displaystyle\end{eqnarray}$$

where $d$ is an unknown function representing the leading-order nonlinear effect in the dispersion relation, which depends on $\unicode[STIX]{x1D70F}$ and the body forcing $f$ . Following the derivation of the linear solution, we further express $A_{g}(\unicode[STIX]{x1D70F},k)$ as

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle A_{g}(\unicode[STIX]{x1D70F},k)=-(k-k_{1g})(k-k_{2g})\quad (k>0), & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle A_{g}(\unicode[STIX]{x1D70F},k)=-(k-k_{3})(k-k_{4})\quad (k<0), & \displaystyle\end{eqnarray}$$

where $k_{1g}$ and $k_{2g}$ are given by the roots of the nonlinear dispersion relation of $A_{g}(\unicode[STIX]{x1D70F},k)=0$ ,

(3.13a,b ) $$\begin{eqnarray}\displaystyle k_{1g}={\textstyle \frac{1}{2}}(1-2\unicode[STIX]{x1D70F}+\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}),\quad k_{2g}={\textstyle \frac{1}{2}}(1-2\unicode[STIX]{x1D70F}-\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}), & & \displaystyle\end{eqnarray}$$

which are the wavenumbers of the two resonant wave modes as $\unicode[STIX]{x1D6FF}\rightarrow 0$ . Here, $k_{3}$ and $k_{4}$ are the wavenumbers of the non-resonance waves, given by

(3.14a,b ) $$\begin{eqnarray}\displaystyle k_{3}={\textstyle \frac{1}{2}}(-1-2\unicode[STIX]{x1D70F}+\sqrt{1+4\unicode[STIX]{x1D70F}}),\quad k_{4}={\textstyle \frac{1}{2}}(-1-2\unicode[STIX]{x1D70F}-\sqrt{1+4\unicode[STIX]{x1D70F}}). & & \displaystyle\end{eqnarray}$$

On substituting $\tilde{h}(k)$ in (3.9) into (3.1) and splitting the integral with $k$ into two parts, we have

(3.15) $$\begin{eqnarray}\displaystyle h(x,z) & = & \displaystyle -\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}\int _{0}^{\infty }\left(\frac{1}{k-k_{1g}}-\frac{1}{k-k_{2g}}\right)\tilde{P}(k)\text{e}^{-\text{i}kx}\text{e}^{|k|z}\,\text{d}k\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{0}\frac{\tilde{P}(k)}{A_{g}(\unicode[STIX]{x1D70F},k)}\text{e}^{-\text{i}kx}\text{e}^{|k|z}\,\text{d}k.\end{eqnarray}$$

Here, both $k_{1g}$ and $k_{2g}$ have positive real parts. In the limit of linear solution (with $d=0$ ), $k_{1g}$ and $k_{2g}$ reduce to $k_{1,2}=(1-2\unicode[STIX]{x1D70F}\pm \unicode[STIX]{x1D6FF})/2$ (e.g. Liu & Yue Reference Liu and Yue1993). At the critical frequency ( $\unicode[STIX]{x1D6FF}=0$ ), $k_{1g,2g}=k_{c}\pm d^{1/2}$ , where the critical wavenumber $k_{c}=1/4$ .

By applying contour integration for the first integral in (3.15) by the use of the indented contours shown in figure 1 and with the imposition of the radiation condition, we obtain

(3.16) $$\begin{eqnarray}\displaystyle h(x,z) & = & \displaystyle -\frac{\text{i}\sqrt{2\unicode[STIX]{x03C0}}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}\tilde{P}(k_{2g})\text{e}^{-\text{i}k_{2g}x+|k_{2g}|z}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\text{i}}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{0}^{\infty }\frac{\tilde{P}(-\text{i}k)}{A(\unicode[STIX]{x1D70F},-\text{i}k)+d}\text{e}^{-kx+|k|z}\,\text{d}k+\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{0}\frac{\tilde{P}(k)}{A(\unicode[STIX]{x1D70F},k)+d}\text{e}^{-\text{i}kx+|k|z}\,\text{d}k\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

for $x>0$ and

(3.17) $$\begin{eqnarray}\displaystyle h(x,z) & = & \displaystyle -\frac{\text{i}\sqrt{2\unicode[STIX]{x03C0}}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}\tilde{P}(k_{1g})\text{e}^{-\text{i}k_{1g}x+|k_{1g}|z}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\text{i}}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{0}^{\infty }\frac{\tilde{P}(\text{i}k)}{A(\unicode[STIX]{x1D70F},\text{i}k)+d}\text{e}^{kx+|k|z}\,\text{d}k+\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{0}\frac{\tilde{P}(k)}{A(\unicode[STIX]{x1D70F},k)+d}\text{e}^{-\text{i}kx+|k|z}\,\text{d}k\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

for $x<0$ . Near the critical frequency, the solution of $h(x,z)$ is dominated by the first term in (3.16) (or (3.17)) for $x>0$ (or $x<0$ ), which is $O(\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FF}^{-1})$ , while the other terms in (3.16) (or (3.17)) are $O(\unicode[STIX]{x1D716})$ in the limit of $d=0$ .

Figure 1.

Sketch of the indented contours and locations of the poles at $k=k_{1g}$ and $k_{2g}$ in the complex $k$ -plane for integrals in (3.15) for (a) $x>0$ and (b) $x<0$ .

To find the unknown function $d$ , we substitute $\tilde{h}(k)$ in (3.9) into (3.7) and carry out the double integration via contour integration. Upon retaining the leading-order (in $\unicode[STIX]{x1D6FF}$ ) contributions due to the residues of the poles of the integrand in the complex $\unicode[STIX]{x1D708}$ - and $\unicode[STIX]{x1D706}$ -planes, we have

(3.18) $$\begin{eqnarray}\displaystyle A(\unicode[STIX]{x1D70F},k)\frac{\tilde{P}(k)}{A_{g}(\unicode[STIX]{x1D70F},k)}-\unicode[STIX]{x03C0}^{2}\frac{a(-k_{c},k_{c},k_{c})\tilde{P}^{\star }(-k_{c})\tilde{P}(k_{c})\tilde{P}(k_{c})}{(\unicode[STIX]{x1D6FF}^{2}+4d)A_{g}(\unicode[STIX]{x1D70F},k)}=\tilde{P}(k), & & \displaystyle\end{eqnarray}$$

where $\tilde{P}^{\star }$ is the Fourier transform of $P^{\ast }(x)$ (cf. (3.3)). From the definition of $\tilde{P}^{\star }$ (cf. (3.3)) and (3.5), it follows that $\tilde{P}^{\star }(-k_{c})=\tilde{P^{\ast }}(k_{c})$ . From (3.8), we have $a(-k_{c},k_{c},k_{c})=-2k_{c}^{4}/\unicode[STIX]{x03C0}$ . By substituting $a(-k_{c},k_{c},k_{c})$ into (3.18) and letting $k=k_{c}$ , we obtain

(3.19) $$\begin{eqnarray}\displaystyle A(\unicode[STIX]{x1D70F},k)+\frac{2\unicode[STIX]{x03C0}k_{c}^{4}\tilde{P}^{\star }(-k_{c})\tilde{P}(k_{c})}{(\unicode[STIX]{x1D6FF}^{2}+4d)}=A_{g}(\unicode[STIX]{x1D70F},k), & & \displaystyle\end{eqnarray}$$

from which we then obtain the equation to determine $d$ in the neighbourhood of the critical frequency,

(3.20) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D6FF}^{2}+4d)d=2\unicode[STIX]{x03C0}k_{c}^{4}|\tilde{P}(k_{c})|^{2}. & & \displaystyle\end{eqnarray}$$

3.2 Solution for a two-dimensional point source

For a two-dimensional point source located at $(0,-H)$ with strength characterized by a small dimensionless parameter $\unicode[STIX]{x1D716}$ , ${\mathcal{X}}$ is given by

(3.21) $$\begin{eqnarray}\displaystyle {\mathcal{X}}(x,z)=\frac{\unicode[STIX]{x1D716}}{\unicode[STIX]{x03C0}}\ln \left(\frac{r}{r_{1}}\right), & & \displaystyle\end{eqnarray}$$

where $r=\sqrt{x^{2}+(z+H)^{2}}$ and $r_{1}=\sqrt{x^{2}+(z-H)^{2}}$ . From (2.11), $P(x)$ is determined as

(3.22) $$\begin{eqnarray}\displaystyle P(x)=-{\mathcal{X}}_{z}(x,0)=\frac{2\unicode[STIX]{x1D716}H}{\unicode[STIX]{x03C0}(x^{2}+H^{2})}. & & \displaystyle\end{eqnarray}$$

Then, we obtain $\tilde{P}(k)$ from (3.5) using (3.22),

(3.23) $$\begin{eqnarray}\displaystyle \tilde{P}(k)=2\unicode[STIX]{x1D716}\frac{\text{i}k}{\sqrt{2\unicode[STIX]{x03C0}}}\text{e}^{-|k|H}. & & \displaystyle\end{eqnarray}$$

From (3.20), $d$ is determined as

(3.24) $$\begin{eqnarray}\displaystyle d=\pm \frac{\unicode[STIX]{x1D716}}{2^{6}}\text{e}^{-k_{c}H}. & & \displaystyle\end{eqnarray}$$

The positive solution of $d$ is consistent with Dagan & Miloh (Reference Dagan and Miloh1982) while the negative solution of $d$ is chosen to be omitted in Dagan & Miloh (Reference Dagan and Miloh1982). Substitution of $d$ in (3.24) into (3.16) and (3.17) leads to $h(x,z)=O(\unicode[STIX]{x1D716}^{1/2})$ in the neighbourhood of the critical frequency (as $\unicode[STIX]{x1D6FF}\rightarrow 0$ ).

3.3 Solution for a two-dimensional body

In the case of a body, ${\mathcal{X}}$ is unknown a priori and it is coupled with $h$ by the body boundary condition (2.9). The potential ${\mathcal{X}}$ needs to be solved together with $h$ by satisfying the free-surface boundary condition (2.8) and the body boundary condition (2.9).

For a two-dimensional submerged body, we express ${\mathcal{X}}(x,z)$ in the form

(3.25) $$\begin{eqnarray}\displaystyle {\mathcal{X}}(x,z)=\int _{\bar{S}_{B}}\unicode[STIX]{x1D70E}(s^{\prime })\ln (r/r_{1})\,\text{d}s^{\prime }, & & \displaystyle\end{eqnarray}$$

where

(3.26a,b ) $$\begin{eqnarray}\displaystyle r=\sqrt{(x-x^{\prime })^{2}+(z-z^{\prime })^{2}},\quad r_{1}=\sqrt{(x-x^{\prime })^{2}+(z+z^{\prime })^{2}} & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D70E}(s^{\prime })$ represents the unknown strength of the source distributed on the body surface. Here, we use $s=(x,z)$ to represent the field point where the potential is calculated while we use $s^{\prime }=(x^{\prime },z^{\prime })$ to represent the location of the source. Substitution of (3.25) into (2.11) yields

(3.27) $$\begin{eqnarray}\displaystyle P(x)=-2\int _{\bar{S}_{B}}\unicode[STIX]{x1D70E}(s^{\prime })\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\ln r\right)_{z=0}\,\text{d}s^{\prime }. & & \displaystyle\end{eqnarray}$$

By using the identity

(3.28) $$\begin{eqnarray}\displaystyle \ln r=\int _{0}^{\infty }\frac{\text{e}^{-k}-\text{e}^{-k|z-z^{\prime }|}\cos k(x-x^{\prime })}{k}\,\text{d}k, & & \displaystyle\end{eqnarray}$$

we can rewrite (3.27) as

(3.29) $$\begin{eqnarray}\displaystyle P(x)=-\int _{-\infty }^{\infty }\text{e}^{-\text{i}kx}\left\{\int _{\bar{S}_{B}}\text{e}^{|k|z^{\prime }}\text{e}^{\text{i}kx^{\prime }}\unicode[STIX]{x1D70E}(s^{\prime })\,\text{d}s^{\prime }\right\}\text{d}k, & & \displaystyle\end{eqnarray}$$

which corresponds to

(3.30) $$\begin{eqnarray}\displaystyle \tilde{P}(k)=-\sqrt{2\unicode[STIX]{x03C0}}\int _{\bar{S}_{B}}\text{e}^{|k|z^{\prime }}\text{e}^{\text{i}kx^{\prime }}\unicode[STIX]{x1D70E}(s^{\prime })\,\text{d}s^{\prime }. & & \displaystyle\end{eqnarray}$$

On substituting $h(x,z)$ in (3.16) and (3.17) as well as ${\mathcal{X}}(x,z)$ in (3.25) into the body boundary condition (2.9), we have

(3.31) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}(s)+\bar{{\mathcal{X}}}_{n}(s)-\frac{\text{i}\sqrt{2\unicode[STIX]{x03C0}}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}(-\text{i}k_{2g}n_{x}+|k_{2g}|n_{z})\tilde{P}(k_{2g})\text{e}^{-\text{i}k_{2g}x+|k_{2g}|z}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{\text{i}}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{0}^{\infty }\left\{\frac{(-kn_{x}+|k|n_{z})\tilde{P}(-\text{i}k)}{A(\unicode[STIX]{x1D70F},-\text{i}k)+d}\text{e}^{-kx+|k|z}\right\}\text{d}k\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{0}\left\{\frac{(-\text{i}kn_{x}+|k|n_{z})\tilde{P}(k)}{A(\unicode[STIX]{x1D70F},k)+d}\text{e}^{-\text{i}kx+|k|z}\right\}\text{d}k=f(s),\quad s(x>0)\in \bar{S}_{B}\end{eqnarray}$$

and

(3.32) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}(s)+\bar{{\mathcal{X}}}_{n}(s)-\frac{\text{i}\sqrt{2\unicode[STIX]{x03C0}}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}(-\text{i}k_{1g}n_{x}+|k_{1g}|n_{z})\tilde{P}(k_{1g})\text{e}^{-\text{i}k_{1g}x+|k_{1g}|z}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\text{i}}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{0}^{\infty }\left\{\frac{(kn_{x}+|k|n_{z})\tilde{P}(\text{i}k)}{A(\unicode[STIX]{x1D70F},\text{i}k)+d}\text{e}^{kx+|k|z}\right\}\,\text{d}k\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\int _{-\infty }^{0}\left\{\frac{(-\text{i}kn_{x}+|k|n_{z})\tilde{P}(k)}{A(\unicode[STIX]{x1D70F},k)+d}\text{e}^{-\text{i}kx+|k|z}\right\}\text{d}k=f(s),\quad s(x<0)\in \bar{S}_{B},\end{eqnarray}$$

where $\bar{{\mathcal{X}}}_{n}(s)$ is defined as

(3.33) $$\begin{eqnarray}\displaystyle \bar{{\mathcal{X}}}_{n}(s)\equiv \unicode[STIX]{x2A0D}_{\bar{S}_{B}}\unicode[STIX]{x1D70E}(s^{\prime })\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}[\ln (r/r_{1})]\,\text{d}s^{\prime }, & & \displaystyle\end{eqnarray}$$

in which the Cauchy principal-value integral is indicated.

To illustrate the characteristic features of the solution near the critical frequency, we here consider the case of a deeply submerged body with centre submergence $H$ . In this case, the two integrals on the left-hand sides of (3.31) and (3.32) behave like ${\sim}O(\text{e}^{-2kH})$ , and they are much smaller than the forcing term $f(s)$ , so that they can be neglected in (3.31) and (3.32). Near the critical frequency, $k_{1g,2g}=k_{c}\pm d^{1/2}\approx k_{c}$ in weak nonlinearity. Using these facts, (3.31) and (3.32) reduce to the same equation,

(3.34) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}(s)+\bar{{\mathcal{X}}}_{n}(s)-\frac{\text{i}\sqrt{2\unicode[STIX]{x03C0}}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}(-\text{i}k_{c}n_{x}+k_{c}n_{z})\tilde{P}(k_{c})\text{e}^{-\text{i}k_{c}x+k_{c}z}=f(s),\quad s\in S_{B}. & & \displaystyle\end{eqnarray}$$

By multiplying each side of (3.34) by the term $\text{e}^{\text{i}k_{c}x+|k_{c}|z}$ and integrating the equation over the body surface $\bar{S}_{B}$ , we obtain

(3.35) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x03C0}\tilde{P}(k_{c})-\text{i}2\unicode[STIX]{x03C0}\frac{\tilde{P}(k_{c})\unicode[STIX]{x1D6E4}k_{c}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}=\sqrt{2\unicode[STIX]{x03C0}}{\mathcal{F}}(k_{c}), & & \displaystyle\end{eqnarray}$$

where

(3.36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}(k_{c})=\int _{\bar{S}_{B}}(-\text{i}n_{x}+n_{z})\text{e}^{2k_{c}z}\,\text{d}s=2k_{c}\int \int _{B}\text{e}^{2k_{c}z}\,\text{d}S & & \displaystyle\end{eqnarray}$$

where $B$ is the mean body section and

(3.37) $$\begin{eqnarray}\displaystyle {\mathcal{F}}(k_{c})=\int _{\bar{S}_{B}}[f(x,z)-\bar{{\mathcal{X}}}_{n}(s)]\text{e}^{\text{i}k_{c}x+k_{c}z}\,\text{d}s. & & \displaystyle\end{eqnarray}$$

From (3.35), we obtain $\tilde{P}(k_{c})$ in terms of $d$ and ${\mathcal{F}}(k_{c})$ ,

(3.38) $$\begin{eqnarray}\displaystyle \tilde{P}(k_{c})=-\!\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{1/2}\frac{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}+\text{i}2k_{c}\unicode[STIX]{x1D6E4}}{\mathcal{F}}(k_{c}). & & \displaystyle\end{eqnarray}$$

Substitution of $\tilde{P}(k_{c})$ in (3.38) into (3.34) yields a new integral equation for $\unicode[STIX]{x1D70E}$ ,

(3.39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}(s)+\bar{{\mathcal{X}}}_{n}(s)+\frac{2k_{c}(n_{x}+\text{i}n_{z}){\mathcal{F}}(k_{c})}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}+\text{i}2k_{c}\unicode[STIX]{x1D6E4}}\text{e}^{-\text{i}k_{c}x+k_{c}z}=f(s),\quad s\in S_{B}. & & \displaystyle\end{eqnarray}$$

It follows from (3.39) that for $\unicode[STIX]{x1D6E4}\neq 0$ , $\unicode[STIX]{x1D70E}=O(\unicode[STIX]{x1D716})$ (since $\unicode[STIX]{x1D70E}(s)\sim f(s)$ ) and is bounded as $\unicode[STIX]{x1D6FF}\rightarrow 0$ .

Substitution of (3.38) into (3.20) gives the equation to determine $d$ ,

(3.40) $$\begin{eqnarray}\displaystyle 4d^{2}+(\unicode[STIX]{x1D6FF}^{2}+4k_{c}^{2}\unicode[STIX]{x1D6E4}^{2})d-4k_{c}^{4}|{\mathcal{F}}(k_{c})|^{2}=0. & & \displaystyle\end{eqnarray}$$

Two solutions for $d$ are obtained as a function of the body forcing $f$ and body geometry parameter $\unicode[STIX]{x1D6E4}$ ,

(3.41) $$\begin{eqnarray}\displaystyle d_{1,2}=-\frac{1}{2}\left(\frac{\unicode[STIX]{x1D6FF}^{2}}{4}+k_{c}^{2}\unicode[STIX]{x1D6E4}^{2}\right)\pm \frac{1}{2}\left[\left(\frac{\unicode[STIX]{x1D6FF}^{2}}{4}+k_{c}^{2}\unicode[STIX]{x1D6E4}^{2}\right)^{2}+4k_{c}^{4}|{\mathcal{F}}(k_{c})|^{2}\right]^{1/2}, & & \displaystyle\end{eqnarray}$$

with $d_{1}$ (or $d_{2}$ ) corresponding to the plus (or minus) sign in (3.41). Once $d$ is determined, we have $\tilde{P}(k_{c})$ from (3.38), $h(x,z)$ from (3.16) and (3.17), and ${\mathcal{X}}(x,z)$ from (3.25), with $\unicode[STIX]{x1D70E}(s)$ from (3.39).

We see from (3.41) that $d_{1}$ is always positive and $d_{1}=O(\unicode[STIX]{x1D716})$ . The wavenumbers of the two resonant wave modes are $k_{1g,2g}=k_{c}\pm (1/2)(\unicode[STIX]{x1D6FF}^{2}+4d_{1})^{1/2}$ , indicating that $k_{1g}>k_{c}$ and $k_{2g}<k_{c}$ . The solution $d_{2}$ is negative definite, and $d_{2}=O(1)$ for $\unicode[STIX]{x1D6E4}\neq 0$ . The associated wavenumbers $k_{1g,2g}=k_{c}\pm (1/2)(\unicode[STIX]{x1D6FF}^{2}+4d_{2})^{1/2}$ become complex as $\unicode[STIX]{x1D6FF}\rightarrow 0$ . In this case, $\text{Im}\,k_{1g,2g}=\pm (1/2)\text{Im}(\unicode[STIX]{x1D6FF}^{2}+4d_{2})^{1/2}$ , which is equal to $\pm |d_{2}|^{1/2}$ at $\unicode[STIX]{x1D6FF}=0$ , representing the spatial damping coefficient of the resonant waves. This solution shows that the two resonant wave modes become evanescent waves near the critical frequency. For $\unicode[STIX]{x1D6E4}=0$ , the solution reduces to $d_{1,2}=O(\unicode[STIX]{x1D716})$ , consistent with the solution for a point source/dipole (Dagan & Miloh Reference Dagan and Miloh1982).

The dominant part of the nonlinear solution of $h(x,z)$ associated with the resonant wave modes can be expressed as

(3.42) $$\begin{eqnarray}\displaystyle h_{NL}(x,z)=\frac{\text{i}2{\mathcal{F}}(k_{c})}{\sqrt{\unicode[STIX]{x1D6FF}^{2}+4d}+\text{i}2k_{c}\unicode[STIX]{x1D6E4}}\text{e}^{-\text{i}k_{1g,2g}x+|k_{1g,2g}|z} & & \displaystyle\end{eqnarray}$$

for $x<0$ and $x>0$ . By setting $d=0$ , equation (3.42) reduces to the linear solution,

(3.43) $$\begin{eqnarray}\displaystyle h_{L}(x,z)=\frac{\text{i}2{\mathcal{F}}(k_{c})}{\unicode[STIX]{x1D6FF}+\text{i}2k_{c}\unicode[STIX]{x1D6E4}}\text{e}^{-\text{i}k_{c}x+k_{c}z}, & & \displaystyle\end{eqnarray}$$

which is consistent with that of Liu & Yue (Reference Liu and Yue1993). From (3.42) and (3.43), we have $h_{NL}=O(\unicode[STIX]{x1D716})$ and $h_{L}=O(\unicode[STIX]{x1D716})$ as $\unicode[STIX]{x1D6FF}\rightarrow 0$ if $\unicode[STIX]{x1D6E4}\neq 0$ . The difference between the amplitudes of $h_{NL}$ and $h_{L}$ characterizes the order of magnitude of the nonlinear correction due to the cubic self-interaction effect,

(3.44) $$\begin{eqnarray}\displaystyle h_{COR1}(x,z)=|h_{NL}|-|h_{L}|\sim \frac{4{\mathcal{F}}(k_{c})\text{e}^{k_{c}z}d_{1}}{\unicode[STIX]{x1D6FF}^{2}+4k_{c}^{2}\unicode[STIX]{x1D6E4}^{2}} & & \displaystyle\end{eqnarray}$$

for $d=d_{1}$ and

(3.45) $$\begin{eqnarray}\displaystyle h_{COR2}(x,z)=|h_{NL}|-|h_{L}|\sim \frac{2{\mathcal{F}}(k_{c})\text{e}^{k_{c}z}(\text{e}^{\pm \unicode[STIX]{x1D708}x}-1)}{\unicode[STIX]{x1D6FF}^{2}+4k_{c}^{2}\unicode[STIX]{x1D6E4}^{2}} & & \displaystyle\end{eqnarray}$$

for $x<0$ and $x>0$ with $d=d_{2}$ , where the spatial damping factor $\unicode[STIX]{x1D708}=\text{Im}(\unicode[STIX]{x1D6FF}^{2}+4d_{2})^{1/2}/2$ . For $\unicode[STIX]{x1D6E4}\neq 0$ , $h_{COR1}=O(\unicode[STIX]{x1D716}^{2})$ while $h_{COR2}=O(\unicode[STIX]{x1D716})$ as $\unicode[STIX]{x1D6FF}\rightarrow 0$ . The nonlinear correction $h_{COR2}$ , associated with $d=d_{2}$ , dominates over $h_{COR1}$ associated with $d=d_{1}$ . Importantly, $h_{COR2}$ is of the same order of magnitude as the (leading-order) linear solution $h_{L}$ . This shows that although the linear solution $h_{L}$ is bounded (for $\unicode[STIX]{x1D6E4}\neq 0$ ), $h_{L}$ is a partial contribution to the leading-order solution. It is necessary to include the cubic nonlinear self-interaction effects in order to obtain the complete leading-order solution near the critical frequency.

If $\unicode[STIX]{x1D6E4}(k_{c})=0$ , the solution reduces to that for a point source with $h_{L}=O(\unicode[STIX]{x1D716}/\unicode[STIX]{x1D6FF})$ and $h_{NL}=O(\unicode[STIX]{x1D716}^{1/2})$ as $\unicode[STIX]{x1D6FF}\rightarrow 0$ , which is consistent with Dagan & Miloh (Reference Dagan and Miloh1982).

We note that $\unicode[STIX]{x1D6E4}$ is always positive definite for a submerged body with non-zero cross-section area (cf. (3.36)). As the cross-section area shrinks to zero, the body becomes like a point and $\unicode[STIX]{x1D6E4}=0$ .

3.4 Extension to a three-dimensional submerged body

The foregoing analysis can be generalized to three dimensions. We define the two-dimensional Fourier transform of $h(x,y,z=0)$ by

(3.46) $$\begin{eqnarray}\displaystyle \tilde{h}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }h(x,y,0)\text{e}^{\text{i}(\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}x+\unicode[STIX]{x1D70C}\sin \unicode[STIX]{x1D703}y)}\,\text{d}x\,\text{d}y. & & \displaystyle\end{eqnarray}$$

By making use of the two-dimensional Fourier transform of (2.8), introducing an unknown small parameter $d$ to the dispersion relation for the nonlinear effects and following the same procedure as in § 3.1, we obtain the dominant part of the solution of $h(x,y,z)$ (for given a free surface forcing $P(x,y)$ ) due to the two resonant wave modes,

(3.47) $$\begin{eqnarray}\displaystyle h(x,y,z)\sim \ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}\tilde{P}(k_{c},0)\text{e}^{-\text{i}k_{1g,2g}x+|k_{1g,2g}z|} & & \displaystyle\end{eqnarray}$$

as $\unicode[STIX]{x1D6FF}\rightarrow 0$ for $x<0$ and $x>0$ , where $\tilde{P}(k_{c},0)$ is the two-dimensional Fourier transform of $P(x,y)$ and $d$ is determined from the equation

(3.48) $$\begin{eqnarray}\displaystyle d\sim [\ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}]^{2}|\tilde{P}(k_{c},0)|^{2}. & & \displaystyle\end{eqnarray}$$

In the presence of a submerged three-dimensional body, we follow the same procedure as in § 3.3 to obtain

(3.49) $$\begin{eqnarray}\displaystyle \tilde{P}(k_{c},0)\sim \frac{{\mathcal{F}}(k_{c})}{\unicode[STIX]{x03C0}+\text{i}k_{c}\unicode[STIX]{x1D6E4}\ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}(k_{c})$ and ${\mathcal{F}}(k_{c})$ are given by (3.36) and (3.37) but with $\bar{S}_{B}$ representing the mean surface of the wetted three-dimensional body surface. Upon substituting (3.49) into (3.47) and (3.48), we obtain the dominant part of the solution for $h(x,y,z)$ associated with the two resonant wave modes,

(3.50) $$\begin{eqnarray}\displaystyle h(x,y,z)\sim \frac{\ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}}{\unicode[STIX]{x03C0}+\text{i}k_{c}\unicode[STIX]{x1D6E4}\ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}}{\mathcal{F}}(k_{c})\text{e}^{-\text{i}k_{1g,2g}x+|k_{1g,2g}z|}, & & \displaystyle\end{eqnarray}$$

and the equation to determine $d$ ,

(3.51) $$\begin{eqnarray}\displaystyle d\sim \frac{[\ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}]^{2}}{\unicode[STIX]{x03C0}^{2}+k_{c}^{2}\unicode[STIX]{x1D6E4}^{2}[\ln (\unicode[STIX]{x1D6FF}^{2}+d)^{1/2}]^{2}}|{\mathcal{F}}(k_{c})|^{2}. & & \displaystyle\end{eqnarray}$$

Two complex-valued roots of $d$ generally exist in (3.51). We see from (3.50) that if $\unicode[STIX]{x1D6E4}\neq 0$ , $h_{NL}=O(\unicode[STIX]{x1D716})$ and $h_{L}=O(\unicode[STIX]{x1D716})$ , which follows from (3.50) by setting $d=0$ . The nonlinear correction, due to dominant cubic self-interactions, is $h_{COR}=|h_{NL}|-|h_{L}|=O(\unicode[STIX]{x1D716})$ , since the complex solutions of $d$ give spatial damping effects on the resonant wave modes, which is the same behaviour as in two dimensions. If $\unicode[STIX]{x1D6E4}=0$ , $h_{NL}=O(\unicode[STIX]{x1D716}\ln \unicode[STIX]{x1D716})$ while $h_{L}=O(\unicode[STIX]{x1D716}\ln \unicode[STIX]{x1D6FF})$ as $\unicode[STIX]{x1D6FF}\rightarrow 0$ , consistent with the point-source solution by Dagan & Miloh (Reference Dagan and Miloh1982).

4.1 Analytic solution

For a submerged circular cylinder with a radius of $R$ and a centre submergence of $H$ , the geometry parameter $\unicode[STIX]{x1D6E4}(k_{c})=2\unicode[STIX]{x03C0}R\text{e}^{-2k_{c}H}I_{1}(2k_{c}R)$ , where $I_{1}$ is the modified Bessel function of the first kind and $\bar{{\mathcal{X}}}_{n}(s)=0$ (Liu & Yue Reference Liu and Yue1993). The function ${\mathcal{F}}(k)$ is then given by

(4.1) $$\begin{eqnarray}\displaystyle {\mathcal{F}}(k)=\int _{\bar{S}_{B}}f(x,z)\text{e}^{\text{i}kx+kz}\,\text{d}s. & & \displaystyle\end{eqnarray}$$

For convenience in the calculation of ${\mathcal{F}}(k)$ , we replace $(x,z)$ by the cylindrical coordinates $(r,\unicode[STIX]{x1D703})$ . For the radiation problem, we have on the body surface

(4.2) $$\begin{eqnarray}\displaystyle f(R,\unicode[STIX]{x1D703})=\unicode[STIX]{x1D709}_{x}\left(\text{i}\unicode[STIX]{x1D70F}\cos \unicode[STIX]{x1D703}+\frac{2\cos 2\unicode[STIX]{x1D703}}{R}\right)+\unicode[STIX]{x1D709}_{z}\left(\text{i}\unicode[STIX]{x1D70F}\sin \unicode[STIX]{x1D703}+\frac{2\sin 2\unicode[STIX]{x1D703}}{R}\right), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D709}_{x}$ and $\unicode[STIX]{x1D709}_{z}$ are respectively the amplitudes of the surge and heave oscillations. Upon expanding the exponential function in (4.1) in a Taylor series, the integration over $\unicode[STIX]{x1D703}$ can be readily carried out to yield

(4.3) $$\begin{eqnarray}\displaystyle {\mathcal{F}}(k_{c})=2\unicode[STIX]{x03C0}k_{c}^{2}R^{2}\text{e}^{-k_{c}H}(-\unicode[STIX]{x1D709}_{x}+\text{i}\unicode[STIX]{x1D709}_{z}). & & \displaystyle\end{eqnarray}$$

Substitution of ${\mathcal{F}}(k_{c})$ in (4.3) into (3.41) gives the solution of $d$ ,

(4.4) $$\begin{eqnarray}\displaystyle d_{2} & = & \displaystyle -\frac{1}{2}\left[\frac{\unicode[STIX]{x1D6FF}^{2}}{4}+4\unicode[STIX]{x03C0}^{2}k_{c}^{2}R^{2}\text{e}^{-4k_{c}H}I_{1}^{2}(2k_{c}R)\right]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{2}\left\{\left[\frac{\unicode[STIX]{x1D6FF}^{2}}{4}+4\unicode[STIX]{x03C0}^{2}k_{c}^{2}R^{2}\text{e}^{-4k_{c}H}I_{1}^{2}(2k_{c}R)\right]^{2}+16k_{c}^{8}\unicode[STIX]{x03C0}^{2}R^{4}\text{e}^{-2k_{c}H}(\unicode[STIX]{x1D709}_{x}^{2}+\unicode[STIX]{x1D709}_{z}^{2})\right\}^{1/2}.\end{eqnarray}$$

The amplitudes of the resonant waves far downstream ( $x<0$ ) and upstream ( $x>0$ ) of the body decay with the distance from the body like

(4.5) $$\begin{eqnarray}\displaystyle A_{1,2}(x)\sim \text{e}^{\unicode[STIX]{x1D708}_{1,2}x} & & \displaystyle\end{eqnarray}$$

for $x<0$ and $x>0$ , with the spatial damping coefficient given by

(4.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}_{1,2}=\pm {\textstyle \frac{1}{2}}\text{Im}(\unicode[STIX]{x1D6FF}^{2}+4d_{2})^{1/2} & & \displaystyle\end{eqnarray}$$

in the neighbourhood of the critical frequency, with $\unicode[STIX]{x1D6FF}\ll 1$ .

4.2 Comparison with direct numerical simulation

To verify the analytical result, we solve the original nonlinear initial boundary-value problem for the submerged circular cylinder in a uniform current by numerical simulation. The comparison of the numerical solution with the analytic prediction provides a verification of the theory.

The simulation is directly performed in the time domain, based on the use of a highly efficient high-order spectral (HOS) method for nonlinear wave–body interactions developed by Liu et al. (Reference Liu, Dommermuth and Yue1992). The time simulation of the nonlinear wave–body interaction problem up to a specified high order $M$ in wave steepness consists of two main steps. Beginning from initial values for the free-surface elevation $\unicode[STIX]{x1D702}(x,t)$ and free-surface potential $\unicode[STIX]{x1D6F7}^{S}(x,t)$ , at each successive time step, we (i) solve the nonlinear boundary-value problem to order $M$ based on the perturbation method to obtain the vertical velocity at the free surface and (ii) integrate the evolution equations, which are the nonlinear kinematic and dynamic free-surface boundary conditions, forward for $\unicode[STIX]{x1D702}(x,t+\unicode[STIX]{x0394}t)$ and $\unicode[STIX]{x1D6F7}^{S}(x,t+\unicode[STIX]{x0394}t)$ , where $\unicode[STIX]{x0394}t$ is the time step; the process is then repeated. By use of the Taylor series expansion of $\unicode[STIX]{x1D6F7}^{S}(x,t)$ (on the instantaneous free surface) about the mean free surface, the nonlinear boundary-value problem at each time is expanded into a sequence of linearized boundary-value problems for the perturbed potentials, $\unicode[STIX]{x1D6F7}^{(m)}(x,z,t)$ , $m=1,\ldots ,M$ , which are successively solved by the pseudo-spectral method with the use of $N_{F}$ Fourier spectral modes on the mean free surface and $N_{B}$ Fourier source modes on the body. With a long-time simulation, we obtain (limited-cycle) steady-state solutions for the surface wave elevation near the body and radiation forces on the body including nonlinear free-surface effects (as in physical laboratory experiments). The harmonic components of the wave elevation and forces on the body are obtained from the steady-state time history via Fourier transform. For more details on the HOS method for nonlinear wave–body simulations, see Liu et al. (Reference Liu, Dommermuth and Yue1992).

We remark that the $N_{F}$ spectral wave modes and $N_{B}$ body-source modes are the basis functions used to efficiently represent the solution of each $\unicode[STIX]{x1D6F7}^{(m)}(x,z,t)$ , $m=1,\ldots ,M$ , at this particular time. The amplitudes of these spectral modes are obtained by invoking the boundary conditions that are not automatically satisfied by each spectral mode. For the nonlinear problem, they are not explicitly related to any specific physical wave or body components in general, and their time dependences are unknown a priori and may not have any physical meaning.

As demonstrated by Liu et al. (Reference Liu, Dommermuth and Yue1992), by using fast-transform techniques, the operational count per time step is (approximately) linearly proportional to $M$ and $N_{F}$ (typically $N_{F}\gg N_{B}$ ). For a submerged body, the method obtains the exponential convergence of the solution with respect to $M$ and $N_{F}$ (for moderately steep waves), and $N_{B}$ . In the present simulation, we use the length of the periodic computational domain $L_{x}=40\unicode[STIX]{x1D706}_{c}$ , $N_{F}=2048$ , $N_{B}=128$ and time step $\unicode[STIX]{x0394}t=T_{c}/128$ , where $\unicode[STIX]{x1D706}_{c}=2\unicode[STIX]{x03C0}/k_{c}$ and the critical wave period $T_{c}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{c}$ . The values of the numerical parameters are determined through convergence tests, with the numerical results unchanged by further varying their values. Simulations with different values of order $M$ are performed to access and quantify the nonlinear free-surface effects. As initial conditions, the free-surface elevation $\unicode[STIX]{x1D702}(x,0)$ and the potential $\unicode[STIX]{x1D6F7}^{S}(x,0)$ are prescribed to be zero. To remove possible reflection from the periodic boundaries, we follow Liu et al. (Reference Liu, Dommermuth and Yue1992) to smoothly truncate the free-surface elevation and velocity potential at longitudinal ends of the computational domain after each time step of the simulation.

Figure 2.

Representative instantaneous free-surface shape of the steady-state wave profile above the submerged circular cylinder (with its centre at $x=0$ ) obtained by the time-domain nonlinear simulation with $M=3$ (——). The analytical solutions of the envelope of decaying wave amplitude of resonant waves downstream ( $x<0$ ) and upstream ( $x>0$ ) of the body, $A_{1,2}(x)\sim \text{e}^{\pm 0.022x}$ , from (4.6) (– – –) are also plotted for comparison. (Here, $H/R=6$ , $F_{r}=0.75$ , $\unicode[STIX]{x1D70F}=0.25$ , $\unicode[STIX]{x1D709}_{x}/R=0.05$ .)

Figure 2 displays the nonlinear simulation result for the steady-state wave profile above the submerged circular cylinder at the critical frequency $\unicode[STIX]{x1D70F}=1/4$ . The nonlinear simulation result (with order $M=3$ ) shows that the resonant wave amplitudes downstream ( $x<0$ ) and upstream ( $x>0$ ) of the body, $A_{1,2}(x)$ , decay with the distance from the body, as predicted by the theory. The analytic solution of the decay rate for the upstream wave amplitude $A_{2}(x)$ compares well with the numerical solution. Downstream, in addition to the resonant wave whose amplitude $A_{1}(x)$ decays with $|x|$ , two non-resonant propagating waves (with wavenumbers $k_{3}$ and $k_{4}$ ) are also present. Unlike upstream of the body, it is difficult to separate the spatially varying amplitude of the resonant wave $A_{1}(x)$ for a clean comparison with the theoretical prediction downstream of the body.

Figure 3.

Comparison of the damping factor $\unicode[STIX]{x1D708}_{2}$ of the resonant wave upstream of the body ( $A_{2}(x)$ ) between the analytic solution from (4.6) with $H/R=4$ (——) and 6 (– – –), and nonlinear simulations ( $M=3$ ) with $H/R=4$ (▴) and 6 (▪) as a function of the surge motion amplitude $\unicode[STIX]{x1D709}_{x}/R$ . (Here, $\unicode[STIX]{x1D70F}=0.25$ and $F_{r}=0.75$ .)

Figure 4.

Comparison of the damping factor $\unicode[STIX]{x1D708}_{2}$ of the resonant wave upstream of the body ( $A_{2}(x)$ ) between the analytic solution from (4.6) at $\unicode[STIX]{x1D70F}=0.245$ (– – –) and 0.25 (——), and nonlinear simulations ( $M=3$ ) at $\unicode[STIX]{x1D70F}=0.245$ (▪) and 0.25 (▴) as a function of body submergence $H/R$ . (Here, $\unicode[STIX]{x1D709}_{x}/R=0.05$ and $F_{r}=0.75$ .)

Figure 3 compares the damping factor of the resonant wave upstream of the body ( $A_{2}(x)$ ) between the analytic solution (4.6) and nonlinear simulations ( $M=3$ ) at the critical frequency $\unicode[STIX]{x1D70F}=1/4$ for different body motion amplitudes. At deep submergence ( $H/R=6$ ), the analytical solution matches well with the nonlinear numerical solution $(M=3)$ . Both the analytic solution and the numerical simulation show that the damping factor $\unicode[STIX]{x1D708}$ has a weak dependence on the body motion amplitude, as indicated in (4.6). In particular, $\unicode[STIX]{x1D708}$ approaches a non-zero asymptote as the body motion amplitude decreases to zero. As the submergence decreases to $H/R=4$ , both the analytic solution and the numerical simulation show the same trend of $\unicode[STIX]{x1D708}$ variation with body motion amplitude, but a quantitative discrepancy exists between them since the analytic solution employs the assumption of deep submergence.

We compare the theoretical solution of $\unicode[STIX]{x1D708}_{2}$ with that from nonlinear numerical simulation ( $M=3$ ) for different body submergences at fixed body motion amplitude. As shown in figure 4, the agreement between the theoretical solution and the nonlinear simulation is excellent at deep submergence for $\unicode[STIX]{x1D70F}=0.245$ and 0.25. The discrepancy between them increases at shallower submergence, as expected, since the theoretical solution is derived based on the assumption of deep submergence. Additionally, both the theoretical prediction and the numerical simulation show that $\unicode[STIX]{x1D708}_{2}$ increases with decreasing body submergence. This is consistent with the fact that $d_{2}\sim \text{e}^{-kH}$ , as shown in (4.4).

The force on the body can be obtained by integration of the pressure over the body in the simulation. The harmonic components can be extracted from the steady-state time history via Fourier transform. The (conventional) added mass and damping coefficients of the body can be determined from the complex amplitude of the first-harmonic force. For example, for surge motion with oscillation $\unicode[STIX]{x1D709}(t)=\unicode[STIX]{x1D709}_{x}\sin (\unicode[STIX]{x1D714}t)$ , the added mass ( $A_{xx}$ and $A_{zx}$ ) and damping ( $B_{xx}$ and $B_{zx}$ ) coefficients of the body can be calculated by

(4.7) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}cc@{}}A_{xx}=2\,\text{Re}\,F_{x}/\unicode[STIX]{x1D709}_{x}\unicode[STIX]{x1D714}^{2}, & A_{zx}=2\,\text{Re}\,F_{z}/\unicode[STIX]{x1D709}_{x}\unicode[STIX]{x1D714}^{2},\\ B_{xx}=-2\,\text{Im}\,F_{x}/\unicode[STIX]{x1D709}_{x}\unicode[STIX]{x1D714}, & B_{zx}=-2\,\text{Im}\,F_{z}/\unicode[STIX]{x1D709}_{x}\unicode[STIX]{x1D714},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $F_{x}$ and $F_{z}$ are the complex amplitudes of the first-harmonic force on the body in the $x$ - and $z$ -directions respectively.

Figure 5.

Time histories of the real part ( $\text{Re}\,F_{x}$ ) and imaginary part ( $\text{Im}\,F_{x}$ ) of the complex amplitude of the first-harmonic radiation force in the horizontal direction on the submerged circular cylinder at $\unicode[STIX]{x1D70F}=0.25$ as a function of simulation time ( $t/T_{c}$ ) obtained by numerical simulations with $M=1$ (●), 2 (▴) and 3 (▪). The results for $\text{Im}\,F_{x}$ with $M=2$ and 3 overlap each other graphically. (Here, $H/R=6$ , $F_{r}=0.75$ , $\unicode[STIX]{x1D709}_{x}/R=0.05$ .)

Figure 6.

The imaginary part (a) and real part (b) of the radiation force in the horizontal direction on the cylinder at $\unicode[STIX]{x1D70F}=0.25$ as a function of the surge motion amplitude ( $\unicode[STIX]{x1D709}_{x}/R$ ) obtained by nonlinear simulations with $M=1$ (●), 2 (▴) and 3 (▪). The results with $M=2$ and 3 overlap each other graphically. (Here, $H/R=6$ and $F_{r}=0.75$ .)

Figure 5 plots the sample time history of the complex amplitude of the first-harmonic horizontal force on the cylinder at the critical frequency $\unicode[STIX]{x1D70F}=1/4$ with a fixed surge motion amplitude, obtained by numerical simulations with $M=1$ , 2 and 3. The results in the figure show that the nonlinear solution (with $M=2$ , 3) approaches the steady state much faster than the linear solution (with $M=1$ ) since the inclusion of nonlinear effects removes the wave resonance. Figure 6 plots the simulation results for the radiation force on the cylinder at the critical frequency $\unicode[STIX]{x1D70F}=1/4$ as a function of the body motion amplitude. The simulation results are converged with $M$ since the difference between the results with $M=2$ and 3 is negligibly small. The results with $M=1$ , 2 and 3, normalized by the body motion amplitude, are all nearly constant and independent of the body motion amplitude. This indicates that the linear solution, nonlinear solution and nonlinear correction due to the cubic self-interactions of resonant waves are all of first order in the wave steepness (or body motion amplitude). This confirms the theoretical prediction derived in § 3.3.

Figure 7.

The damping coefficient (a) and added mass coefficient (b) of the submerged circular cylinder in the neighbourhood of $\unicode[STIX]{x1D70F}=0.25$ obtained by nonlinear simulations with $M=1$ (●), 2 (▴) and 3 (▪). The results with $M=2$ and 3 overlap each other graphically. (Here, $H/R=6$ , $F_{r}=0.75$ and $\unicode[STIX]{x1D709}_{x}/R=0.05$ .)

Figure 7 shows the numerical simulation results with $M$ up to 3 for the wave damping and added mass coefficients of the body in the neighbourhood of $\unicode[STIX]{x1D70F}=1/4$ . For both wave damping and added mass, the linear solution $(M=1)$ is finite but varies sharply near $\unicode[STIX]{x1D70F}=1/4$ , while the nonlinear solution ( $M=3$ ) varies smoothly, with the peak value near $\unicode[STIX]{x1D70F}=1/4$ significantly reduced. Such important nonlinear corrections are seen to exist only in a small region near $\unicode[STIX]{x1D70F}=1/4$ . For $\unicode[STIX]{x1D70F}$ away from $1/4$ , the nonlinear effects become insignificant. The results in figures 6 and 7 indicate that the nonlinear corrections due to the third-order self-interactions should be accounted for in order to obtain the correct prediction of the seakeeping solution near the critical frequency.

The analysis in § 3.3 shows that the behaviour of the solution for a body reduces to that of a point source/dipole in the neighbourhood of $\unicode[STIX]{x1D70F}=1/4$ when the geometry parameter $\unicode[STIX]{x1D6E4}=0$ . For a submerged body, the value of $\unicode[STIX]{x1D6E4}$ vanishes exponentially with increasing submergence. We thus use the numerical simulation result for a deeply submerged circular cylinder ( $H/R=20$ ) to check the validity of the theory of Dagan & Miloh (Reference Dagan and Miloh1982) for a point source/dipole, which has not been achieved in the literature. The linear simulation (with $M=1$ ) at $\unicode[STIX]{x1D70F}=1/4$ does not reach the steady state, with the solution growing with time indefinitely, indicating that the linear solution is singular in this case. The nonlinear simulation (with $M=3$ ) reaches the steady state rapidly. Figure 8 plots the (steady-state) normalized amplitude of the radiated wave ( $A_{2}/\unicode[STIX]{x1D709}_{x}$ ) upstream of the cylinder ( $x/R=80$ ) at $\unicode[STIX]{x1D70F}=1/4$ as a function of body motion amplitude ( $\unicode[STIX]{x1D709}_{x}/R$ ). The nonlinear simulations also show that $A_{2}$ does not decay with the distance from the body in this case. The result in figure 8 shows that $A_{2}/\unicode[STIX]{x1D709}_{x}$ behaves as $(\unicode[STIX]{x1D709}_{x}/R)^{-1/2}$ , which indicates the nonlinear solution of $A_{2}\sim (\unicode[STIX]{x1D709}_{x}/R)^{1/2}=O(\unicode[STIX]{x1D716}^{1/2})$ . These numerical results for a deeply submerged body (with $\unicode[STIX]{x1D6E4}\sim 0$ ) confirm the theory of Dagan & Miloh (Reference Dagan and Miloh1982) for a point source/dipole.

Figure 8.

The amplitude of the radiated wave upstream of the deeply submerged circular cylinder, $A_{2}/\unicode[STIX]{x1D709}_{x}$ , as a function of the surge motion amplitude, $\unicode[STIX]{x1D709}_{x}/R$ . The plotted results are from nonlinear simulations with $M=3$ (▪) and a fitted dependence $(\unicode[STIX]{x1D709}_{x}/R)^{-1/2}$ (– – –). (Here, $\unicode[STIX]{x1D70F}=0.25$ , $H/R=20$ and $F_{r}=0.75$ .)