2.1. Shared governing equations

For each of the three flow configurations considered in this paper, we assume the fluid is incompressible and inviscid and that the flow is irrotational. We choose to ignore surface tension to specifically observe the effects of a low Froude number on gravity waves (as opposed to gravity–capillary waves). In a real world context, our observations would apply to large cargo ships, near-surface submarines manoeuvring slowly or slow-moving currents past floating platforms, for example, where Froude numbers can be small while surface tension is negligible. On the hand, our analysis may not be directly applicable for slowly moving small-scale experimental hulls or a duck, for example, where surface tension is important.

We denote the location of the unknown free surface by $z=\zeta (x,y)$. The velocity potential $\phi (x,y,z)$ therefore satisfies Laplace's equation throughout the flow field

(2.1)\begin{equation} \nabla^2\phi=0,\quad\text{for }z<\zeta(x,y).\end{equation}

In dimensional units we suppose the disturbance (submerged source, submerged doublet or pressure) is moving with speed $U$ and is associated with a representative length scale $L$. We set $g$ to be acceleration due to gravity. Thus, by fixing our frame of reference to move with the disturbance, the dimensionless kinematic and dynamic conditions on the free surface are

(2.2)\begin{gather} \phi_x\zeta_x+\phi_y\zeta_y=\phi_z,\quad\text{on }z=\zeta(x,y), \end{gather}

(2.3)\begin{gather}(\phi_x^2+\phi_y^2+\phi_z^2)+\frac{2\zeta}{F^2}+\delta p(x,y)=1,\quad\text{on }z=\zeta(x,y), \end{gather}

where the Froude number is defined by

(2.4)\begin{equation} F=\frac{U}{\sqrt{gL}}. \end{equation}

The term $\delta p(x,y)$ is an applied pressure on the surface (in addition to atmospheric pressure), as described below. The appropriate far-field conditions for our dimensionless problems are

(2.5)\begin{gather} \phi\sim x,\quad\zeta\rightarrow 0\quad \text{as }x\rightarrow-\infty, \end{gather}

(2.6)\begin{gather}\phi\sim x,\quad\text{as }z\rightarrow-\infty. \end{gather}

Equation (2.5) enforces the radiation condition that surface gravity waves do not propagate ahead of the disturbance, while (2.6) simply ensures that the flow approaches a uniform stream in the infinitely deep limit.

2.2. Flow past a submerged point source

The first flow configuration involves a point source of dimensional strength $m$, submerged a dimensional distance $L$ from the (undisturbed) surface, moving steadily through a fluid in a horizontal direction with dimensional speed $U$. The motion of the source produces a steady ship wake pattern on the free surface. We fix our frame of reference to move with the source, so that our governing equations are (2.1)–(2.6), except that we do not have an additional pressure on the surface, so we set $\delta =0$ in (2.3). Further, we have the additional condition

(2.7)\begin{equation} \phi\sim x-\frac{\epsilon}{4{\rm \pi}\sqrt{x^2+y^2+(z+1)^2}}\quad\text{as }(x,y,z)\rightarrow(0,0,-1),\end{equation}

where $\epsilon =m/(UL^2)$, which simply ensures that the velocity potential has the appropriate singular behaviour at the source.

The linearised version of this problem ($\epsilon \ll 1$; see § 3.1) is equivalent to flow past a submerged semi-infinite Rankine body with a rounded nose (like a cigar) which approaches a cylinder of radius $\sqrt {\epsilon /{\rm \pi} }$ in the far field (Batchelor Reference Batchelor1967). Nonlinearity acts to distort the shape of the submerged body, and for highly nonlinear flows ($\epsilon \gg 1$), the analogy with a submerged Rankine body no longer holds (Pethiyagoda et al. Reference Pethiyagoda, McCue and Moroney2014a). In what follows, the numerical solutions of the fully nonlinear equations (2.1)–(2.6) (with $\delta =0$) and (2.7) are computed using the boundary integral method outlined in Pethiyagoda et al. (Reference Pethiyagoda, McCue, Moroney and Back2014b), which is based on the algorithm outlined in Forbes (Reference Forbes1989).

2.3. Flow past a submerged point doublet

For our second flow configuration, we replace the point source singularity in § 2.2 with a point doublet of dimensional strength $\kappa$. Thus, the governing equations are again (2.1)–(2.6) with $\delta =0$, but instead of (2.7) we have

(2.8)\begin{equation} \phi\sim x+\frac{\mu x}{4{\rm \pi}(x^2+y^2+(z+1)^2)^{3/2}},\quad\text{as }(x,y,z)\rightarrow(0,0,-1),\end{equation}

where $\mu =\kappa /(UL^3)$ is dimensionless strength of the doublet.

In this case, the linear version of the problem ($\mu \ll 1$; see § 3.2) is equivalent to flow past a submerged solid sphere of radius $(\mu /2{\rm \pi} )^{1/3}$ (Lamb Reference Lamb1916). Nonlinear solutions for moderate values of $\mu$ have a very similar interpretation, except that the spherical body is distorted and indeed the surface is no longer closed (Pethiyagoda et al. Reference Pethiyagoda, McCue and Moroney2014a). All of our numerical solutions to the nonlinear problem (2.1)–(2.6), (2.8) ($\delta =0$) are computed using the scheme outlined in Pethiyagoda et al. (Reference Pethiyagoda, McCue and Moroney2014a).

2.4. Flow past a pressure distribution

Our third configuration involves studying the wake that forms behind a steadily moving pressure distribution applied to the free surface. The type of pressure we are focussing on is localised and characterised by some pressure scale $P$ and horizontal length scale $L$. By fixing our frame of reference to move with the pressure, the steady problem is to solve (2.1)–(2.6), where $\delta =P/\rho U^2$ is the dimensionless pressure strength and $\rho$ is fluid density.

In our formulation, the localised pressure distribution could take any form (Miao & Liu Reference Miao and Liu2015), provided it decays to zero as $x^2+y^2\rightarrow \infty$. For the most part we use the simple Gaussian

(2.9)\begin{equation} p(x,y)=\exp({-{\rm \pi}^2(x^2+y^2)}),\end{equation}

which has been used extensively in the past as a very simple model for generating ship wakes (Darmon et al. Reference Darmon, Benzaquen and Raphaël2014; Ellingsen Reference Ellingsen2014; Pethiyagoda et al. Reference Pethiyagoda, McCue and Moroney2015; Li & Ellingsen Reference Li and Ellingsen2016; Pethiyagoda, McCue & Moroney Reference Pethiyagoda, McCue and Moroney2017; Smeltzer & Ellingsen Reference Smeltzer and Ellingsen2017; Pethiyagoda et al. Reference Pethiyagoda, Moroney, Macfarlane, Binns and McCue2018). In the present form, (2.1)–(2.6), (2.9) together describe a nonlinear free-surface flow problem which depends on the Froude number $F$ and the pressure strength $\delta$. In §§ 4–6 we present numerical results for this problem which are calculated using the numerical scheme outlined in Pethiyagoda et al. (Reference Pethiyagoda, McCue and Moroney2017), which is based on earlier work by Părău & Vanden-Broeck (Reference Părău and Vanden-Broeck2002) and Părău, Vanden-Broeck & Cooker (Reference Părău, Vanden-Broeck and Cooker2005).

2.5. Measuring the apparent wake angle

We consider two different methods of measuring the apparent wake angle of a ship wake, $\theta _{app}$. Method I was used by Darmon et al. (Reference Darmon, Benzaquen and Raphaël2014), where a wave crest is isolated downstream and then the wave elevation is plotted against the polar angle as measured from the origin to the centreline ($y=0$). We choose an arbitrary percentage of the maximum height and mark the polar angle that gives this height as the apparent wake angle $\theta _{app}$. For this paper we will use 5 %, 10 % and 20 % in our measurements ($\alpha =0.05$, $0.1$ and $0.2$). We use Method I for low Froude numbers where the wake is dominated by transverse waves.

The second method of measuring the apparent wake angle $\theta _{app}$, Method II, given in Pethiyagoda et al. (Reference Pethiyagoda, McCue and Moroney2014a), is performed by isolating each transverse wavelength of the wake and marking the highest peak. We then fit a line to the highest peaks using a linear least squares algorithm. The apparent wake angle $\theta _{app}$ is then given by the angle between the fitted line and the centreline. We use Method II for moderate values of the Froude number, for which wakes are made up of both transverse and divergent waves.

3.1. Submerged point source ($\epsilon \rightarrow 0$ with $F$ fixed)

Taking the limit as $\epsilon \rightarrow 0$ while fixing $F$ has the effect of turning off the source while keeping its speed constant. We can linearise the problem about this limit to give a modified system of governing equations, namely Laplace's equation

(3.1)\begin{equation} \nabla^2\phi=0,\quad\text{for }z<0,\end{equation}

together with the kinematic condition

(3.2)\begin{equation} \zeta_x=\phi_z \quad \text{on} \ z=0,\end{equation}

and the dynamic condition

(3.3)\begin{equation} \phi_x-1+\frac{\zeta}{F^2}=0 \quad \text{on} \ z=0.\end{equation}

The remaining governing equations (2.5)–(2.7) remain unchanged. Note that the kinematic and dynamic boundary conditions have been projected onto the plane $z=0$, which is a (well-known) key feature of this linearisation.

We will use the exact solution to (3.1)–(3.3), (2.5)–(2.7), given by Peters (Reference Peters1949), namely

(3.4)\begin{align} \zeta(x,y)&={-}\frac{\epsilon F^2 \,\mathrm{sgn}(x)}{{\rm \pi}^2}\int_{0}^{{\rm \pi}/{2}}\cos\psi \int_{0}^{\infty}\frac{k \,\mathrm{e}^{{-}k|x|}\cos(k y \sin \psi)g(k,\psi)}{F^4 k^2+\cos^2\psi}\, \mbox{d}k\,\mbox{d}\psi \nonumber\\ &\quad +\frac{\epsilon H(x)}{\rm \pi}\int_{-\infty}^{\infty}\xi \,\mathrm{e}^{{-}F^2\xi^2} \cos(x\xi)\cos(y\xi \lambda)\,\mbox{d}\lambda, \end{align}

where

(3.5)\begin{gather} g(k,\psi)=F^2 k \sin(k\cos \psi)+\cos\psi\cos(k\cos\psi), \end{gather}

(3.6)\begin{gather}\xi(\lambda)=\sqrt{\lambda^2+1}/F^2, \end{gather}

and $H(x)$ is the Heaviside function. In this form, the double integral provides the near-field component that decays rapidly away from the origin, while the single integral provides the wave train component. In figure 1(a–c) we show the wave pattern from the single integral for three different Froude numbers. For $F=0.1$ and $0.2$, we can see the wake is completely dominated by transverse waves, while for $F=1$ we see a mix of both transverse and divergent waves. Note that the reason we only show the wave pattern in figure 1 from the single integral is that, for small Froude number $F$, the near-field contribution from the double integral is much larger than the wave amplitude, and so when we plot the full surface the waves are obscured by the near-field disturbance.

Figure 1. (a–c) Plan view of the wave pattern for linear flow past a submerged source for Froude numbers $F=0.1$, $0.2$ and $1$. Only the single integral in (3.4) is used for the computation. (d–f) Equivalent images for flow past a submerged doublet, this time computed using the single integral term of (3.16). In (a,b) and (d,e), the dashed line indicates the wake angle defined by 20 % of the wave height. In (c,f), the dashed line is Kelvin's angle $\theta _{{wedge}}$. (g) A plot of wave height against angle ($\theta$) for flow past a point source (solid curves) and a doublet (dashed curves) given by (3.15) and (3.19), respectively, for the Froude numbers $F=0.2$ (blue, dashed orange), $F=0.1$ (green, dashed violet) and $F=0.05$ (red, dashed blue).

We perform a stationary phase approximation on the exact linear solution (3.4) which results in the approximate wave profile given in polar coordinates ($r=\sqrt {x^2+y^2}$, $\theta =\tan ^{-1}(y/x)$)

(3.7)\begin{equation} \zeta(r,\theta) \sim a_1(r,\theta)\cos\left(r g(\lambda_1(\theta),\theta)+\frac{\rm \pi}{4}\right) +a_2(r,\theta)\cos\left(r g(\lambda_2(\theta),\theta)-\frac{\rm \pi}{4}\right) \end{equation}

as $r\rightarrow \infty$ for $|\theta |<\arcsin (1/3)$, where

(3.8a,b)\begin{gather} \lambda_1(\theta)=\frac{-1+\sqrt{1-8\tan^2\theta}}{4\tan\theta},\quad \lambda_2(\theta)=\frac{-1-\sqrt{1-8\tan^2\theta}}{4\tan\theta}, \end{gather}

(3.9a,b)\begin{gather} f(\lambda)=\frac{\epsilon\sqrt{\lambda^2+1}}{F^2}\exp\left({-\frac{\lambda^2+1}{F^2}}\right), \quad g(\lambda,\theta)=\frac{\sqrt{\lambda^2+1}}{F^2} (\cos\theta+\lambda\sin\theta), \end{gather}

(3.10a,b)\begin{gather} a_1(r,\theta)=\sqrt{\frac{2}{\rm \pi}}\frac{f(\lambda_1(\theta))}{\sqrt{r|g_{\lambda\lambda} (\lambda_1(\theta),\theta)|}}, \quad a_2(r,\theta)=\sqrt{\frac{2}{\rm \pi}}\frac{f(\lambda_2(\theta))} {\sqrt{r|g_{\lambda\lambda} (\lambda_2(\theta),\theta)|}}, \end{gather}

and $g_{\lambda \lambda }(\lambda ,\theta )$ is the second partial derivative of $g(\lambda ,\theta )$ with respect to $\lambda$. In simple terms, the far-field surface takes the form

(3.11)\begin{equation} \zeta(r,\theta)\sim \frac{\epsilon}{r^{1/2}}A(\theta)\exp({-G(\theta)/F^2}) \quad\mbox{as}\ r\rightarrow\infty, \end{equation}

where $A$ and $G$ are functions of $\theta$ (Pethiyagoda et al. Reference Pethiyagoda, McCue and Moroney2014a). The linear dependence on $\epsilon$ in (3.11) arises because we are considering the linear version of the problem. The $r^{-1/2}$ dependence implies that the wave pattern decays rather slowly as we move further away from the cause of the waves. Finally, the exponential dependence on $-F^{-2}$ reveals that the waves themselves become extremely small in the small-Froude-number limit, which is the limit we are most interested in here.

We now consider (3.7) in more detail. The factors $a_1$ and $a_2$ represent the function envelopes for the transverse and divergent waves, respectively. Since $\lambda _1^2<\lambda _2^2$ for $0\leqslant \theta < \theta _{{wedge}}$ (where remember that $\theta _{{wedge}}$ is Kelvin's angle), it follows that $\exp ({-(\lambda _2^2+1)/{F^2}})\ll \exp ({-(\lambda _1^2+1)/{F^2}})$ for $F\ll 1$. In other words, for small Froude numbers, the transverse waves exponentially dominate the divergent waves and thus, in practice, we will observe the highest parts of the wave pattern along the crestlines of the transverse waves.

The location of the crestlines is given by

(3.12)\begin{equation} \cos\left(r g(\lambda_1(\theta),\theta)+\frac{\rm \pi}{4}\right)=1, \end{equation}

which leads to the relationship

(3.13)\begin{equation} r(\theta)=\frac{(2n-1/4){\rm \pi}}{g(\lambda_1(\theta),\theta)},\end{equation}

where $n$ is some positive integer. Thus, we have a formula for the wave height along the crest

(3.14)\begin{equation} a^s_{crest}(\theta)=\epsilon\sqrt{\frac{2g(\lambda_1(\theta),\theta)}{(2n-1/4){\rm \pi}^2g_{\lambda\lambda} (\lambda_1(\theta),\theta)}}\,f(\lambda_1(\theta)). \end{equation}

Scaling such that $\hat {a}^s_{crest}(0)=1$ gives

(3.15)\begin{align} \hat{a}^s_{crest}(\theta)=(\lambda_1(\theta)^2+1)^{3/2} \sqrt{\frac{\cos\theta+\lambda_1(\theta)\sin\theta}{\cos\theta+(2\lambda_1(\theta)^3+3\lambda_1(\theta)) \sin\theta}}\exp({-{\lambda_1(\theta)^2}/{F^2}}).\end{align}

Plotting this crest height for different Froude numbers (figure 1g), we see that as the Froude number $F$ decreases, the Gaussian-like wave crest narrows, which means that the apparent wake angle decreases as $F$ decreases. We can use (3.15) to determine the wake angle, $\theta _{app}$, by setting $\hat {a}^s_{crest}(\theta _{app})=\alpha$, where $\alpha$ is the chosen fraction of the maximum wave height. For example, in figure 1(a,b) we have plotted dashed lines that indicate the apparent wake angle when $\alpha$ is chosen to be $\alpha =0.2$.

3.2. Submerged point doublet ($\mu \rightarrow 0$ with $F$ fixed)

The linear governing equations for flow past a submerged doublet are given by (3.1)–(3.3), (2.5)–(2.6) and (2.8). They are derived by fixing $F$ in the fully nonlinear equations and taking the limit $\mu \rightarrow 0$, which is equivalent to considering flow past a submerged sphere of vanishing radius while keeping its depth and the far-field speed constant. The exact solution to the linear problem is given by

(3.16)\begin{align} \zeta(x,y)&=\frac{\mu F^2}{{\rm \pi}^2}\int_{0}^{{\rm \pi}/{2}}\cos\psi \int_{0}^{\infty}\frac{k^2 \,\mathrm{e}^{{-}k|x|}\cos(k y \sin \psi)g(k,\psi)}{F^4 k^2+\cos^2\psi}\, \mbox{d}k\,\mbox{d}\psi \nonumber\\ &\quad -\frac{\mu H(x)}{\rm \pi}\int_{-\infty}^{\infty}\xi^2 \,\mathrm{e}^{{-}F^2\xi^2}\sin(x\xi)\cos(y\xi \lambda)\, \mbox{d}\lambda, \end{align}

where $g(k,\psi )$ and $\xi (\lambda )$ are given by (3.5) and (3.6), respectively. Representative wave patterns computed using the single integral term in (3.16) are shown in figure 1(d–f). For the two small Froude numbers $F=0.1$ and $0.2$, the plan view looks very similar to that for flow past a submerged source, with the wake appearing to be made up of transverse waves only, and the waves themselves appearing to be a phase shift when compared to flow past a submerged source. For the moderate Froude number $F=1$, the image is more like a traditional ship wake, with both transverse and divergent waves.

As with the flow due to a submerged source, we can write the stationary phase approximation

(3.17)\begin{equation} \zeta(r,\theta) \sim a_1(r,\theta)\sin\left(r g(\lambda_1(\theta),\theta)+\frac{\rm \pi}{4}\right) +a_2(r,\theta)\sin\left(r g(\lambda_2(\theta),\theta)-\frac{\rm \pi}{4}\right), \end{equation}

where $a_1(r,\theta )$, $a_2(r,\theta )$, $\lambda _1(\theta )$, $\lambda _2(\theta )$ and $g(\lambda ,\theta )$ are given by (3.8a,b), (3.9b)–(3.10b), but this time

(3.18)\begin{equation} f(\lambda)=\mu(\lambda^2+1)\exp({-{(\lambda^2+1)}/{F^2}}).\end{equation}

Following the same argument, the scaled amplitude along the transverse wave crest is given by

(3.19)\begin{equation} \hat{a}^d_{crest}(\theta)=(\lambda_1(\theta)^2+1)^2\sqrt{\frac{\cos\theta+\lambda_1(\theta)\sin\theta}{ \cos\theta+(2\lambda_1(\theta)^3+3\lambda_1(\theta)) \sin\theta}} \exp({-{\lambda_1(\theta)^2}/{F^2}}).\end{equation}

We note that $\hat {a}^d_{crest}(\theta )=\sqrt {(\lambda _1(\theta )^2+1)}\hat {a}^s_{crest}(\theta )$, and $\lambda _1\rightarrow 0$ as $F\rightarrow 0$; therefore, the shape of the transverse wave crest for flow due to a submerged doublet will approach that due to a submerged source in the low-Froude-number limit. This property can be seen in figure 1(g), where the curves are actually different but virtually indistinguishable on this scale.

3.3. Applied pressure distribution ($\delta \rightarrow 0$ with $F$ fixed)

The linear governing equations for flow past a pressure distribution are given by (3.1) and (3.2), (2.5) and (2.6) and the linearised dynamic condition

(3.20)\begin{equation} \phi_x-1+\frac{\zeta}{F^2}+p(x,y)=0 \quad \text{on} \ z=0,\end{equation}

where the applied pressure distribution is (2.9). One form for the exact solution, adapted from Pethiyagoda et al. (Reference Pethiyagoda, McCue and Moroney2017), is

(3.21)\begin{align} \zeta(x,y) &={-}\delta F^2 p(x,y)\nonumber\\ &\quad +\frac{\delta F^2}{2{\rm \pi}^3} \int_{-{\rm \pi}/2}^{{\rm \pi}/2}\,\int_{0}^{\infty}\frac{k^2\exp({-k^2/4{\rm \pi}^2}) \cos(k[|x|\cos\psi+y\sin\psi])}{k-k_0}\,\mathrm{d}k\,\mathrm{d}\psi\nonumber\\ &\quad -\frac{\delta F^2 H(x)}{{\rm \pi}^2}\int_{-\infty}^{\infty}\xi^2 \exp({-F^4\xi^4/4{\rm \pi}^2})\sin(x\xi)\cos(y\xi \lambda)\,\mbox{d}\lambda, \end{align}

where the path of $k$-integration is diverted below the pole $k=1/(F^2\cos ^2\psi )$.

As with the submerged doublet, the stationary phase approximation is given by (3.17), where $a_1(r,\theta )$, $a_2(r,\theta )$, $\lambda _1(\theta )$, $\lambda _2(\theta )$ and $g(\lambda ,\theta )$ are given by (3.8a,b), (3.9b)–(3.10b),

(3.22)\begin{equation} f(\lambda)=\delta F^2(\lambda^2+1)\exp({-{(\lambda^2+1)^2}/{4{\rm \pi}^2F^4}})\end{equation}

and the scaled transverse wave crest is given by

(3.23)\begin{align} \hat{a}^p_{crest}(\theta) &=(\lambda_1(\theta)^2+1)^2\sqrt{\frac{\cos\theta+\lambda_1(\theta)\sin\theta}{\cos\theta+(2\lambda_1(\theta)^3+3\lambda_1(\theta))\sin\theta}}\nonumber\\ &\quad \times \exp({-{((\lambda_1(\theta)^2+1)^2-1)}/{4{\rm \pi}^2F^4}}). \end{align}

As mentioned above, we are using the applied pressure distribution (2.9) to align this work with a number of previous studies. It is interesting to note that the choice

(3.24)\begin{equation} p(x,y)={-}\frac{\partial}{\partial x}\left(\frac{1}{2{\rm \pi}\sqrt{x^2+y^2+1}}\right), \end{equation}

provides a wave pattern that is identical to that produced by a submerged source. Similarly,

(3.25)\begin{equation} p(x,y)={-}\frac{\partial^2}{\partial x^2}\left(\frac{1}{2{\rm \pi}\sqrt{x^2+y^2+1}}\right), \end{equation}

provides a wave pattern that is identical to that produced by a submerged doublet (McCue, Pethiyagoda & Moroney Reference McCue, Pethiyagoda and Moroney2019). Further, if we apply one derivative in (3.25) using the product rule, then one of the resulting two terms corresponds to a pressure distribution considered by Havelock (Reference Havelock1919). This pressure leads to a different near-field wave pattern to that due to flow past a submerged doublet, however the far-field behaviour is the same.

3.4. Apparent wake angle for linear flows

In figure 2 we plot apparent wake angles for linearised flow past a submerged source (blue), submerged doublet (red) or an applied pressure distribution (violet). The results using Method I (see § 2.5), valid for low Froude numbers, are represented by the solid curves for a 20 % cutoff height ($\alpha =0.2$). The open circles, triangles and squares are computed using Method II. The dashed lines show the large-Froude-number approximation for Method II found using the method of stationary phase.

Figure 2. Plots of apparent wake angle against the Froude number. The solid curve represents Method I defined in § 2.5 using a 20 % cutoff height. The low-Froude-number results only consider the transverse waves of the stationary phase approximation. The apparent wake angles calculated from Method II are given by the marks with the large-Froude-number approximation given by the dashed line. The blue, red and yellow colours denote flow past a submerged source, submerged doublet and applied pressure distribution, respectively.

We see from figure 2 that, for low Froude numbers, the apparent wake angle $\theta _{app}$ appears to depend linearly on the Froude number $F$ on this log scale, which suggests a power-law relationship. For flow past a point source we can determine the relationship by first setting the desired amplitude $\hat {a}^s_{crest}=\alpha$, taking logs of both sides and assuming $\theta _{app}\ll 1$ to give

(3.26)\begin{equation} \ln\alpha=\left(\frac{5}{2}-\frac{1}{F^2}\right)\theta^2+{O}(\theta^4).\end{equation}

Rearranging (3.26) for $\theta_{app}$ and taking $F\rightarrow 0$, we find

(3.27)\begin{equation} \theta_{app}\sim \sqrt{\ln(1/\alpha)}\,F\quad\text{as }F\rightarrow 0. \end{equation}

Performing the same procedure to flow past a doublet will give the same scaling (3.27), while for flow past the pressure distribution (2.9) we find

(3.28)\begin{equation} \theta_{app}\sim {\rm \pi}\sqrt{2\ln(1/\alpha)}\,F^2\quad\text{as }F\rightarrow 0. \end{equation}

In a sense, figure 2 combines the results of Darmon et al. (Reference Darmon, Benzaquen and Raphaël2014) and Pethiyagoda et al. (Reference Pethiyagoda, McCue and Moroney2014a) (which are for moderate to high Froude numbers) and fills in the gaps by including our new results for smaller Froude numbers.

We note that the ${O}(F^2)$ scaling for flow past the pressure distribution (2.9) is different from the ${O}(F)$ scaling for flow past a submerged source/doublet, other applied pressure distributions can also give rise to the ${O}(F)$ result, or indeed another scaling entirely. To clarify this issue, in appendix A we show how the wake angle scaling for a general pressure distribution depends heavily on the Fourier transform of the pressure and its derivatives.

3.5. Wake envelope via exponential asymptotics

In our analysis above we have applied the method of stationary phase to determine the far-field wake and, subsequently, the apparent wake angle $\theta _{app}$, for linear flow past a submerged point source, submerged doublet and applied pressure distribution. In this subsection we concentrate on flow past a submerged source or doublet, and apply ideas from exponential asymptotics to derive an envelope function which encloses the low-Froude-number wake in both the near and far field.

We suspect from our stationary phase argument that the wave amplitude for the linear problems is of order $\mathrm {e}^{-\chi /F^2}$ for $F\ll 1$, where $\chi$ is a function of position. A straightforward asymptotic expansion implies that, for small $F$, the surface height can be represented by

(3.29)\begin{equation} \zeta(x,y)=\sum_{n=0}^{N-1} F^{2n} \zeta_n(x,y) + Z_N(x,y), \end{equation}

where $Z_N(x,y)$ is the remainder after $N$ terms. If the summation in (3.29) is taken to infinity, the series is divergent; however, if $N$ is chosen carefully so that the series is truncated optimally (at the least term) (Boyd Reference Boyd1999), then

(3.30)\begin{equation} Z_N \sim B(x,y) \,\mathrm{e}^{-\chi/F^2}\quad\mbox{as}\ F\rightarrow 0.\end{equation}

In other words, in order to capture the wave pattern whose amplitude is smaller than all orders of the traditional asymptotic series in powers of $F^2$, we must analyse the remainder term after the original series is truncated optimally. A leading-order expression for $Z_N$ is called a superasymptotic approximation (Berry Reference Berry1991).

The key issue here is that the singulant function $\chi$ governs how small $Z_N$ will be. In particular, as $\mathrm {Re}(\chi )$ increases, then $Z_N$ decreases exponentially. Thus, we are motivated to plot curves of constant $\mathrm {Re}(\chi )$, with the goal of deriving an effective envelope (or boundary) of the wake for small Froude number. For the geometry of flow past a submerged point source, this problem has been studied by Lustri & Chapman (Reference Lustri and Chapman2013). We will not repeat the details here as they are extensive, except to say that in order to derive a partial differential equation (PDE) for $\chi$, they write out an ansatz for the terms $\zeta _n$ in the limit $n\rightarrow$ and (together with an analogous expression for the late-order terms in a series for $\phi$) substitute this into the governing equations to a coupled system of PDEs. These in turn combine to give

(3.31)\begin{equation} \left(\frac{\partial\chi}{\partial x}\right)^4+\left(\frac{\partial\chi}{\partial x}\right)^2+\left(\frac{\partial\chi}{\partial y}\right)^2=0 ,\end{equation}

subject to the condition that $\chi$ must vanish at the singularity of the analytically continued free surface ($x,y \in \mathbb {C}, z =0$), namely

(3.32)\begin{equation} \chi = 0 \quad \mathrm{on} \ x^2 + y^2 + 1 = 0.\end{equation}

Using Charpit's method, this PDE can be solved to give

(3.33)\begin{equation} \chi={\pm}\frac{s-x}{s(2+s^2)},\end{equation}

where $s$ is one of the four roots to the quartic

(3.34)\begin{equation} (x^2+y^2)s^4+4x s^3+(x^2+4y^2+4)s^2+4x s+(4y^2+4)=0.\end{equation}

We are interested in the curves defined by

(3.35)\begin{equation} \mathrm{Re}(\chi)=\nu, \end{equation}

where $\nu$ is some constant, the positive value of $\chi$ and $s$ such that $\mathrm {Re}(x s)>0$ is chosen. There is no need to distinguish between complex conjugates.

To relate (3.35) to the apparent wake angle $\theta _{app}$ we consider the leading-order asymptotic approximation of $\mathrm {Re}(\chi )$ for $r\rightarrow \infty$ (in polar coordinates), derived in appendix B,

(3.36)\begin{equation} \mathrm{Re}(\chi)=f(\theta)+{O}(r^{{-}1}),\end{equation}

where

(3.37)\begin{align} f=\frac{3\cos\theta\sqrt{9\cos^2\theta-8}-9\cos^2\theta+8} {27\cos^4\theta-9\cos^3\theta\sqrt{9\cos^2\theta-8}-42\cos^2\theta+10\cos\theta\sqrt{9\cos^2\theta-8}+16}.\end{align}

Therefore, for some constant $\nu$, the curve $\mathrm {Re}(\chi )=\nu$ approaches a ray of angle $\theta$ measured from the origin defined as the solution to $f(\theta )=\nu$. In figure 3 we present contour plots, where the contours are given for every 10th percentile of the wave height within a single wavelength, for the free-surface profiles given by the single integral term of (3.4). We also include the curve defined by $\mathrm {Re}(\chi )=\nu$ where for figure 3(a,b) $f(\theta _{app})=\nu$, $\hat {a}_{crest}(\theta _{app})=\alpha$ and $\alpha =0.1$. For figure 3(c,d) the constant $\nu$ is chosen by ensuring the curve passes through the edge of the 10 % ($\alpha =0.1$) contour. We can see that the curve (3.35) closely borders the 10th percentile contour. Note that we have drawn plots like those in figure 3 using different values of $\alpha$ (not shown here) and the agreement is just as good as that shown in figure 3.

Figure 3. A contour plot of the positive heights of the free-surface profile for flow past a submerged source, in which there is a contour for every 10th percentile of the wave height within a single wavelength. The solid line is a line of constant $\mathrm {Re}(\chi )$ as given in (3.33), specifically (a) $\mathrm {Re}(\chi )= 1.0236$, (b) $\mathrm {Re}(\chi )= 1.1026$, (c) $\mathrm {Re}(\chi )= 1.4170$ and (d) $\mathrm {Re}(\chi )= 2.3602$. The dashed line is Kelvin's angle $\theta _{{wedge}}$ measured from the origin.

The image in figure 3(a), which is for a very small Froude number $F=0.1$, focuses on the near-field wave pattern. To appreciate the sense in which the curve $\mathrm {Re}(\chi )=\nu$ bends around and approaches a constant ray in the far field, we show in figure 4 the same image but drawn on a much larger scale. Taking together, both figures 3(a) and 4 demonstrate how well the curve $\mathrm {Re}(\chi )=\nu$ acts to define the envelope of the wave pattern for small Froude numbers.

Figure 4. A contour plot of the positive heights of the free-surface profile for flow past a source with $F=0.1$, The contour is at 10 % of the wave height within a single wavelength. The solid line is given by $\mathrm {Re}(\chi )= 1.0236$. The dashed line is Kelvin's angle $\theta _{{wedge}}$.

We emphasise that the analysis in Lustri & Chapman (Reference Lustri and Chapman2013) was performed for the problem of linearised flow past a submerged point source, using the condition (2.7). However, the analysis used to derive the governing equation for $\chi$ (3.31) and the boundary condition (3.32) does not depend on the behaviour or strength of the singularity at $z=-1$, but only its location. Therefore, the expression for $\chi$ given in (3.33) may be applied in an identical fashion to derive an envelope function for the analogous problem of flow past a submerged doublet (or even a range of other point singularities). We illustrate this behaviour in figure 5, where lines of constant $\mathrm {Re}(\chi )$ are overlaid on the positive heights of the free-surface profile for flow past a submerged doublet with $F = 0.1$ and $0.2$. We see that these curves provide an accurate representation of the envelope of the apparent wave region on the surface.

Figure 5. A contour plot of the positive heights of the free-surface profile for flow past a submerged doublet, in which there is a contour for every 10th percentile of the wave height within a single wavelength. The solid line is a line of constant $\mathrm {Re}(\chi )$ as given in (3.33), specifically (a) $\mathrm {Re}(\chi )= 1.0236$, (b) $\mathrm {Re}(\chi )= 1.1026$. The dashed line is Kelvin's angle $\theta _{{wedge}}$.

3.6. Interference effects at low Froude numbers

When analysing the apparent wake angle for ship wakes with large Froude numbers, Noblesse et al. (Reference Noblesse, He, Zhu, Hong, Zhang, Zhu and Yang2014) and Zhang et al. (Reference Zhang, He, Zhu, Yang, Li, Zhu, Lin and Noblesse2015) have shown that interference between bow and stern waves have the effect of changing the large-Froude-number scaling in this limit. Here we will briefly address the effects of interference on the apparent wake angle in the low-Froude-number regime. In figure 6(a–c) we present surface profiles for flow past a Rankine body constructed by a submerging a point source and point sink of equal strengths, where the sink is downstream a horizontal distance $\ell$ from the source (here the dimensionless depth of the source and sink is unity). In this way, the bow-like waves created by flow past the submerged source may interfere with the stern-like waves due to flow past the submerged sink. The solutions we present in figure 6, generated from the superposition of surface profiles given by (3.4) with $F=0.2$, are shown for separation distances (a) $\ell =1.89$, (b) $\ell =1.97$ and (c) $\ell =2.01$, representing maximum constructive, intermediate and destructive interference, respectively.

Figure 6. (a–c) Plan view of the wave pattern for linear flow past a submerged source and sink, of equal strengths, separated a non-dimensional horizontal distance of (a) $\ell =1.89$ (constructive interference), (b) $\ell =1.97$ (no net interference) and (c) $\ell =2.01$ (destructive interference) for Froude number $F=0.2$. A superposition of the single integral in (3.4) is used for the computation. The dashed line indicates the wake angle defined by $\alpha =0.2$ for each of the point singularities. The surface height along the red curves in (a–c) are plotted in (d) against ray angle $\theta$.

The values of $\ell$ chosen for figure 6 come from noting that the wavelength of a transverse wave is $2{\rm \pi} F^2$; for destructive interference we require that $\ell /2{\rm \pi} F^2$ is an integer, while for constructive interference the same quotient is a half-integer. The intermediate interference profile was chosen such that the quotient, $\ell /2{\rm \pi} F^2$, has a remainder of $\pm 1/6$. Note the associated length-based Froude numbers $F_\ell =F/\sqrt {\ell }$ are (a) $F_\ell =0.1457$, (b) $F_\ell =0.1425$ and (c) $F_\ell =0.1410$. Using colour intensity to denote wave amplitude, it is easy to see from this figure that even though the three geometries are very similar, the effects of wave interference on wave amplitude can be significant. Indeed, the waves in (a) appear to have a much larger amplitude than in (c). On the other hand, the apparent wake angle does not appear to change from panels (a–c).

To better understand this comparison, we have also presented in figure 6(d) the surface elevation along crests marked by the red curves in figure 6(a–c). In this figure we have normalised the amplitudes so that the intermediate wave in (b) has a height of unity at $\theta =0$. It can be readily seen in figure 6(d) that the maximum elevation of a wave crest varies greatly depending on whether there is constructive or destructive interference between the bow-lie and stern-like waves. The different maximum elevations mean Method I will give very different apparent wake angle measurements, depending on what type of interference is occurring. Therefore, using Method I to measure the apparent wake angle will no longer give a monotonically decreasing wake angle when Froude number decreases because the surface profile will cycle through constructive and destructive interference as the Froude number changes.

While the maximum height between the three crests in figure 6(d) vary, all three crest collapse onto the same curve as the ray angle moves away from the centreline. This is because the downstream sink only effects the surface profile within a wedge that is contained wholly with then a wedge that defines the influence of the submerged source. As can be seen in figure 6(a–c), between the dashed lines is a region on the outer edge of the wake pattern for which only the upstream source contributes to the surface profile. Therefore, even though Method I will give different results between a Rankine body and a single submerged source, the wake angle measurement for a source is applicable for a Rankine body provided its length is sufficiently large.

4.1. Submerged point source ($F\rightarrow 0$ with $\epsilon$ fixed)

Taking the limit of $F\rightarrow 0$ with $\epsilon$ fixed is analogous to flow past a submerged Rankine body with a rounded nose (that approaches a cylinder of radius $\sqrt {\epsilon /{\rm \pi} }$ in the far field) while reducing the speed of the base flow to zero. In order to keep the dimensionless parameter $\epsilon$ constant, the dimensional source strength $m$ must also be decreased such that the ratio $m/U$ is constant, thus in this interpretation the source is being turned off as its speed is slowed.

In the asymptotic regime $F\rightarrow 0$ with $\epsilon$ fixed, we can approximate the free-surface profile by performing a simple rigid lid approximation similar to that provided in Forbes & Hocking (Reference Forbes and Hocking1990). To this end, we write

(4.1a,b)\begin{equation} \phi = \phi_0(x,y,z)+{O}(F^2),\quad \zeta = F^2\zeta_1(x,y)+{O}(F^4), \end{equation}

and perturb $\phi$ about $z=0$ to give

(4.2)\begin{align} \phi(x,y,F^2\zeta_1+{O}(F^4))= \phi_0(x,y,0)+F^2\left[\phi_1(x,y,0)+\zeta_1(x,y)\phi_{0z}(x,y,0)\right]+{O}(F^4).\end{align}

We then substitute (4.1a,b) and (4.2) into the kinematic and dynamic conditions to derive appropriate boundary conditions for the rigid lid approximation.

The first-order term of the kinematic condition gives us $\phi _{0z}(x,y,0)=0$. Using the method of images with the condition (2.7), we find

(4.3)\begin{equation} \phi_0(x,y,z)=x-\frac{\epsilon}{4{\rm \pi}}\left(\frac{1}{\sqrt{x^2+y^2+(z+1)^2}}+\frac{1}{\sqrt{x^2+y^2+(z-1)^2}}\right). \end{equation}

Taking the first-order terms of the dynamic condition and substituting (4.3) in yields

(4.4)\begin{equation} \zeta_1(x,y)={-}\frac{\epsilon^2}{8{\rm \pi}^2}\frac{x^2+y^2}{(x^2+y^2+1)^3}-\frac{\epsilon}{2{\rm \pi}}\frac{x}{(x^2+y^2+1)^{3/2}}. \end{equation}

From (4.1a,b), the leading-order approximation is $\zeta =F^2\zeta _1$.

The rigid lid solution (4.4) is not capable of capturing any of the details of the wave pattern behind the submerged body; however, it does approximate the near-field behaviour extremely well, at least for small Froude numbers. To make such a comparison, we have in figure 7 representative free-surface profiles for the example $F=0.3$, $\epsilon =0.1$. In figure 7(a), the numerical solution to the fully nonlinear problem (2.1)–(2.6) (with $\delta =0$) and (2.7) is shown, while in figure 7(b) we have the rigid lid solution (4.1b) with (4.4). On this scale, both surfaces appear to compare very well with each other.

Figure 7. A comparison of three different free-surface profiles for flow past a point source with $F=0.3$ and $\epsilon =0.1$. The different profiles are of; (a) the full nonlinear solution, (b) the linear solution, (c) the rigid lid approximation and (d) the difference between the nonlinear solution and the rigid lid approximation. In order to give a sense of the scale, the maximum surface heights in these images are (a) $6.27\times 10^{-4}$, (b) $5.50\times 10^{-4}$, (c) $6.27\times 10^{-4}$ and (d) $1.00\times 10^{-4}$. The nonlinear solution is presented on a $721\times 241$ mesh with $\Delta x=\Delta y=0.05$ and $x_0=-18$.

Also included in figure 7 is the linear solution (3.4) from § 3.1. For these chosen parameters, the linear solution is also very similar to the full numerical solution, although the match is not as good as the rigid lid solution. The one key difference is that the linear solution is able to capture the small transverse waves, while (as just mentioned) the rigid lid solution is not.

Given the amplitude of the waves downstream for the disturbance is extremely small for this example, it is difficult to extract these waves from the full numerical solution in panel (a) of figure 7. We can, however, subtract the rigid lid solution (4.1b), (4.4) from the numerical solution to isolate these waves, as in figure 7(d). Thus the rigid lid solution turns out to be useful to help analyse the full numerical solution for small Froude numbers.

In order to perform a full asymptotic analysis in the limit $F\rightarrow 0$ with $\epsilon$ fixed, we require techniques in exponential asymptotics that have not yet been fully developed. We discuss this issue in § 7.

4.2. Submerged point doublet ($F\rightarrow 0$ with $\mu$ fixed)

By fixing $\mu$ and taking the limit $F\rightarrow 0$, the configuration tends to a flow of decreasing magnitude past a submerged spherical body (of radius $(\mu /2{\rm \pi} )^{1/3}$). In dimensional terms, fixing $\mu$ while decreasing $F$ is equivalent reducing both $U$ and $\kappa$ while keeping $\kappa /U$ constant. Thus, another interpretation of this limit is the doublet is gradually being turned off as it is being slowed down.

The rigid lid approximation can also be applied to this flow configuration. Omitting the details, we find

(4.5)\begin{align} \zeta(x,y)&=F^2\left[\frac{\mu^2}{8{\rm \pi}^2}\left(-\frac{1}{(x^2+y^2+1)^3}+\frac{6x^2}{(x^2+y^2+1)^4}-\frac{9(x^4+x^2y^2)}{(x^2+y^2+1)^5}\right)\right.\nonumber\\ &\quad\left.+\frac{\mu}{2{\rm \pi}}\left(-\frac{1}{(x^2+y^2+1)^{3/2}}+\frac{3x^2}{(x^2+y^2+1)^{5/2}}\right)\right]. \end{align}

For brevity, we shall not include images of the surface (4.5); however, it turns out that the results are analogous to figure 7 in that the rigid lid solution (4.5) provides a very good approximation for the full nonlinear solution in the near field, but is deficient in the sense that it does not describe the wave train downstream. Further discussion on this asymptotic limit is deferred until § 7.

4.3. Applied pressure distribution ($F\rightarrow 0$ with $\delta$ fixed)

Here we suppose that $F\rightarrow 0$ with $\delta$ fixed, which has the physical interpretation of turning the pressure off as it is being slowed down. The rigid lid approximation for this regime is given by the straightforward expression

(4.6)\begin{equation} \zeta(x,y)={-}F^2\delta p(x,y). \end{equation}

Again, to save space we shall not go into any further details here, except to repeat that the rigid lid approximation does a very good job of predicting the near-field behaviour for small $F$, but does not capture the wave pattern, as the wave amplitudes are exponentially small compared to powers of $F^2$.

5.1. Submerged point source ($F\rightarrow 0$ with $\epsilon F$ fixed)

In the double limit $F\rightarrow 0$, $\epsilon \rightarrow \infty$ with $\epsilon F$ held constant, the solution approaches that for flow due to a stationary submerged point source explored in Forbes & Hocking (Reference Forbes and Hocking1990), Vanden-Broeck & Keller (Reference Vanden-Broeck and Keller1997) and Hocking, Forbes & Stokes (Reference Hocking, Forbes and Stokes2016), governed by the single parameter $\bar {m}=\epsilon F=m/\sqrt {gL^5}$. This base flow is itself a highly non-trivial free boundary problem, which these authors demonstrate has solutions up to some limiting value $\bar {m}$, $\bar {m}_{{max}}$ say. There is some uncertainty about what this limiting value is, given the challenges introduced by the nonlinearities involved. In figure 8(a), we present an axisymmetric free-surface profile for this base case for the single value $\bar {m}=2$. This is a highly nonlinear solution which is close the limiting configuration for $\bar {m}=\bar {m}_{{max}}$.

Figure 8. A free-surface profile for nonlinear flow past (a–c) a point source with $\epsilon F=2$, (d–f) a point doublet with $\mu F=1$ and (g–i) a pressure distribution with $\delta F^2=3$. All solutions are presented on a $721\times 361$ mesh with $\Delta x=\Delta y=0.033$ and $x_0=-12$.

To illustrate the regime with $F\ll 1$, $\bar {m}$ constant, we include two nonlinear free-surface profiles in figure 8(b,c) with $\bar {m}=2$, this time for $F=0.03$ ($\epsilon =200/3$) and $F=0.07$ ($\epsilon =200/7$). These are clearly non-axisymmetric surfaces which appear to be small perturbations of the solution in figure 8(a). They correspond to slowly moving submerged point sources.

Note that in figure 8(a) there is a central peak in the surface elevation at the origin. This peak is accompanied by a stagnation point which must occur on the surface directly above the source due to symmetry. In panels (b,c), for which there is a non-zero flow in the far field, we see that the central peak has moved slightly along the $x$-axis in the negative direction. In these two cases, the height of the peaks is less than $F^2/2$ and so these peaks are not stagnation points.

5.2. Submerged point doublet ($F\rightarrow 0$ with $\mu F$ fixed)

If we take $F\rightarrow 0$ and $\mu \rightarrow \infty$ with $\mu F=\bar {\mu }$ held constant, then the limiting flow is that due to a stationary point doublet. As far as we can tell, the fully nonlinear version of this configuration does not appear to be considered in the literature. We have not undertaken a detailed analysis of the parameter space for this interesting problem, as it is outside the scope of the present study. Instead, we present a representative profile in figure 8(d) for $\bar {\mu }=1$. Then in figure 8(e,f) we show solutions with the same value of $\bar {\mu }=1$, but with $F=0.02$ and $0.04$, respectively. These profiles with $F>0$ are clearly small perturbations of the base configuration with a stationary point doublet. They correspond to a slowly moving submerged point doublet.

5.3. Applied pressure distribution ($F\rightarrow 0$ with $\delta F$ fixed)

In the limit $F\rightarrow 0$, $\delta \rightarrow \infty$ with $\delta F^2=\bar {\delta }$ fixed, the flow approaches that due to a stationary pressure distribution. An example of this trivial configuration is shown in figure 8(g), where $\bar {\delta }=3$. We can perturb this flow by moving the applied pressure distribution very slowly, for example as in figure 8(h,i). These two rather unremarkable solutions are also for $\bar {\delta }=3$, except this time $F=0.06$ and $0.12$, respectively.

6.1. Submerged point source ($\epsilon \rightarrow \epsilon _{{max}}$ with $F$ fixed)

To begin, consider the problem of flow past a submerged point source. By setting $F=0.3$, $\epsilon =0.1$, we can compute a solution which is only weakly nonlinear, as shown in figure 9(a). Now, as $\epsilon$ increases, the effects of nonlinearity are more prominent. A train of waves appears behind the disturbance. If we keep increasing $\epsilon$, then we can attempt to explore the limit $\epsilon \rightarrow \epsilon _{{max}}$ with $F$ fixed, where $\epsilon _{{max}}$ is defined such that the free surface $\zeta (x,y)$ for $\epsilon =\epsilon _{{max}}$ has the property $\max (\zeta )=F^2/2$ (as permitted by Bernoulli's equation (2.3)). This would be the limiting configuration characterised by a stagnation point where $\zeta (x,y)=\max (\zeta )$. With our numerical scheme, we are not able to compute the actual limiting configuration, but we are able to resolve solutions that are close to the limit. For example, with $F=0.3$, the largest value of $\epsilon$ we were able to compute a solution was $\epsilon =8.5$. The corresponding free-surface profile is shown in figure 9(b). As an indication of how nonlinear this solution is, the maximum slope of the surface along the centreline is $32^\circ$, which occurs immediately before the highest peak. It is interesting to note this steepness is slightly higher than that for the steepest possible two-dimensional wave ($30^\circ$) (Schwartz Reference Schwartz1974), known as the Stokes limiting configuration (Stokes Reference Stokes1847).

Figure 9. Free-surface profiles for flow past (a,b) a submerged point source, (c,d) a submerged point doublet and (e,f) an applied pressure distribution for Froude number $F=0.3$. Both weakly nonlinear and (a,c,e) highly nonlinear (b,d,f) solutions shown. The maximum gradient of the surface elevation along the centreline for each case is approximately (a) $0.02^\circ$, (b) $32^\circ$, (c) $0.3^\circ$, (d) $25^\circ$, (e) $5^\circ$ and (f) $35^\circ$. All solutions are presented on a $721\times 241$ mesh with $\Delta x=\Delta y=0.05$ and $x_0=-18$.

Looking closely at figure 9(a,b), we see there is an initial peak in the free-surface profile. This initial peak is the only prominent feature in figure 9(a), for which the strength of the source is small. On the other hand, in figure 9(b), at the highest source strength we could solve for, there is also the initial peak, but this is followed by many further peaks along the wave train. It turns out that the first peak of the wave train is higher than the initial peak in figure 9(b). If we assume the solution will reach its limiting configuration when $\zeta (x,y)=\max (\zeta )$ at one of the peaks, it is of interest to speculate which peak is the cause of the solution breaking down as $\epsilon \rightarrow \epsilon _{{max}}$.

We do this in figure 10(a,b) for two values of the Froude number, where we plot the height of the initial peak and the first peak of the wave train against the strength of the source, $\epsilon$, remembering that the higher the value of $\epsilon$ the more nonlinear the solution is. We also plot the limiting height $F^2/2$ as a horizontal dashed line, remembering that no solutions exist with any value of $\zeta >\max (\zeta )=F^2/2$. For $F=0.2$, in figure 10(a) we see the blue data points representing the initial peak are much higher than the first peak of the wave train, and so it is almost certain that the initial peak is approaching the limiting height as $\epsilon \rightarrow \epsilon _{{max}}$. In contrast, for $F=0.3$ the pattern appears the same for small and moderate values of $\epsilon$, with the initial peak much higher than the first peak of the wave train, but then for a very large value of $\epsilon$ the data intersect. Therefore, for $F=0.3$ it appears that the first peak of the wave train is the cause of the solution reaching the limiting configuration. As far as we can tell, while this type of issue is explored in some detail for steady flows in two dimensions, the more complicated analogues in three dimensions are relatively unexplored in the literature (Pethiyagoda et al. Reference Pethiyagoda, McCue and Moroney2014a; Buttle et al. Reference Buttle, Pethiyagoda, Moroney and McCue2018).

Figure 10. A plot of the height of the initial peak (blue) and the first peak of the wave train (red) against the key parameter ($\epsilon$, $\mu$ or $\delta$) for flow past (a,b) a submerged source, (c,d) a submerged doublet and (e,f) an applied pressure distribution. Plots in (a,c,e) are for $F=0.2$ while those in (b,d,f) are for $F=0.3$. The dashed line is the Stokes limiting height, $\zeta =F^2/2$.

6.2. Submerged point doublet ($\mu \rightarrow \mu _{{max}}$ with $F$ fixed)

For the problem of flow past a submerged point doublet, there are some notable differences when compared to the previous subsection. First, we show a weakly nonlinear solution in figure 9(c) for $\mu =0.1$ and a highly nonlinear solution in figure 9(d) for $\mu =1.58$. While there are some relatively small waves behind the disturbance for $\mu =0.1$, these waves dominate the wake for $\mu =1.58$. The solution presented in figure 9(d) is the most nonlinear solution we could compute in the sense that $\mu =1.58$ is the largest value of the doublet strength for which our numerical scheme would converge. These waves here are very steep, with the maximum slope along the centreline roughly $25^\circ$. A clear qualitative difference between this highly nonlinear solution and the one in figure 9(b) for flow past a point source is that the wave pattern in figure 9(d) is characterised by very steep divergent waves. This is surprising as the Froude number here is very low and a rule of thumb is that ship wave profiles for small Froude number are dominated by transverse waves.

Another observation from figure 9(c,d) is that in both cases we see the initial peak of the free-surface profile is lower than the first peak of the wave train, and in the highly nonlinear case in figure 9(d) the difference in heights is substantial. This means that, ultimately, the solution will reach a limiting configuration as $\mu \rightarrow \mu _{{max}}$ when the height of the first peak of the wave train reaches $\max (\zeta )=F^2/2$. To demonstrate how these two peaks depend on the strength of the doublet, we show data in figure 10(c,d) for two values of the Froude number $F=0.2$ and $0.3$. The qualitative difference between figures 10(c,d) and 10(a,b) is clear.

6.3. Applied pressure distribution ($\delta \rightarrow \delta _{{max}}$ with $F$ fixed)

We have also studied the problem of flow past a pressure distribution by fixing the Froude number $F$ and increasing the strength of the pressure $\delta$. We compare two free-surface profiles in figure 9(e,f). The solution in figure 9(e) is only weakly nonlinear as $\delta$ is rather small, while the solution in figure 9(f) is highly nonlinear with maximum slope along the centreline of roughly $35^\circ$. As with figure 9(d), we see in figure 9(f) that the wave pattern is dominated by the divergent waves, which is surprising since the Froude number is small and we may expect the divergent waves to be much smaller than the transverse waves in this regime. The strong nonlinearity here is obviously causing this counterintuitive behaviour.

It is interesting to also track the height of the first peak of the wave train as $\delta$ increases. For $F=0.2$ we see in figure 10(e) that this first peak appears to be increasing with $\delta$ in a way that suggests the solution will reach a limiting configuration when this maximum height reaches $\max (\zeta )$ as $\delta \rightarrow \delta _{{max}}$. On the other hand, for $F=0.3$ in figure 10(f), the trend is not as convincing, and so it is not clear whether a further increase in $\delta$ will lead to the solution ultimately reaching a limiting configuration.