2.1 Problem definition

We consider the irrotational flow of an inviscid, incompressible fluid bounded above by a free surface and below by a channel bottom with arbitrary topography. Far upstream, there is a uniform stream flowing with constant velocity $U$ and depth $H$ . Cartesian coordinates are set up such that the flow is predominantly in the positive $x$ -direction and gravity acts in the negative $z$ -direction. Due to the bottom topography, a wave pattern resembling a ship’s wake appears on the free surface. A schematic of the problem set-up is shown in figure 1. We are interested in the steady-state problem that arises from the long-time limit of this set-up.

Figure 1.

Diagram of the flow problem. The $x$ -axis is aligned with the direction of the flow far upstream. The fluid is bounded above by the free surface $z=\unicode[STIX]{x1D701}(x,y)$ and bounded below by the topography $z=\unicode[STIX]{x1D6FD}(x,y)$ , which, in this schematic, is an isolated bump.

The problem is non-dimensionalised by scaling all lengths by the upstream depth, $H$ , and all velocities by $U$ . By denoting the free surface by $z=\unicode[STIX]{x1D701}(x,y)$ and the bottom boundary by $z=\unicode[STIX]{x1D6FD}(x,y)$ , the dimensionless problem is to solve Laplace’s equation for the velocity potential $\unicode[STIX]{x1D6F7}(x,y,z)$ ,

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6F7}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}^{2}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}^{2}}{\unicode[STIX]{x2202}z^{2}}=0,\quad \unicode[STIX]{x1D6FD}(x,y)<z<\unicode[STIX]{x1D701}(x,y),\end{eqnarray}$$

subject to the kinematic, dynamic and radiation conditions:

(2.2) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6F7}_{x}\unicode[STIX]{x1D701}_{x}+\unicode[STIX]{x1D6F7}_{y}\unicode[STIX]{x1D701}_{y}=\unicode[STIX]{x1D6F7}_{z},\quad & & \displaystyle z=\unicode[STIX]{x1D701}(x,y),\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6F7}_{x}\unicode[STIX]{x1D6FD}_{x}+\unicode[STIX]{x1D6F7}_{y}\unicode[STIX]{x1D6FD}_{y}=\unicode[STIX]{x1D6F7}_{z},\quad & & \displaystyle z=\unicode[STIX]{x1D6FD}(x,y),\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle \displaystyle \frac{1}{2}(\unicode[STIX]{x1D6F7}_{x}^{2}+\unicode[STIX]{x1D6F7}_{y}^{2}+\unicode[STIX]{x1D6F7}_{z}^{2})+\frac{\unicode[STIX]{x1D701}}{F^{2}}=\frac{1}{2},\quad & & \displaystyle z=\unicode[STIX]{x1D701}(x,y),\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle \displaystyle (\unicode[STIX]{x1D6F7}_{x},\unicode[STIX]{x1D6F7}_{y},\unicode[STIX]{x1D6F7}_{z})\rightarrow (1,0,0),\quad & & \displaystyle x\rightarrow -\infty ,\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D701}\rightarrow 0,\quad & & \displaystyle x\rightarrow -\infty ,\end{eqnarray}$$

where $F$ is the depth-based Froude number, $F=U/\sqrt{gH}$ . We assume the dimensionless bottom topography $z=\unicode[STIX]{x1D6FD}(x,y)$ has the property $\unicode[STIX]{x1D6FD}\rightarrow -1$ as $x^{2}+y^{2}\rightarrow \infty$ ; that is, the bottom surface becomes flat away from the origin.

The velocity potential on the free surface $z=\unicode[STIX]{x1D701}(x,y)$ is denoted by $\unicode[STIX]{x1D719}(x,y)=\unicode[STIX]{x1D6F7}(x,y,\unicode[STIX]{x1D701}(x,y))$ ; on the bottom boundary $z=\unicode[STIX]{x1D6FD}(x,y)$ , it is denoted by $\unicode[STIX]{x1D713}(x,y)=\unicode[STIX]{x1D6F7}(x,y,\unicode[STIX]{x1D6FD}(x,y))$ . With this notation, the dynamic condition (2.4) can be rewritten to incorporate the kinematic condition (2.2) as

(2.7a,b ) $$\begin{eqnarray}\displaystyle \frac{1}{2}\frac{(1+\unicode[STIX]{x1D701}_{x}^{2})\unicode[STIX]{x1D719}_{y}^{2}+(1+\unicode[STIX]{x1D701}_{y}^{2})\unicode[STIX]{x1D719}_{x}^{2}-2\unicode[STIX]{x1D701}_{x}\unicode[STIX]{x1D701}_{y}\unicode[STIX]{x1D719}_{x}\unicode[STIX]{x1D719}_{y}}{1+\unicode[STIX]{x1D701}_{x}^{2}+\unicode[STIX]{x1D701}_{y}^{2}}+\frac{\unicode[STIX]{x1D701}}{F^{2}}=\frac{1}{2},\quad z=\unicode[STIX]{x1D701}(x,y). & & \displaystyle\end{eqnarray}$$

We use this form of the boundary condition below.

At this stage it is worth noting that once a solution is computed, the pressure acting on the channel bottom can be calculated by rearranging and evaluating Bernoulli’s equation along the bottom boundary to give

(2.8) $$\begin{eqnarray}p(x,y,\unicode[STIX]{x1D6FD})=\frac{1}{2}-\frac{1}{2}\frac{(1+\unicode[STIX]{x1D6FD}_{x}^{2})\unicode[STIX]{x1D713}_{y}^{2}+(1+\unicode[STIX]{x1D6FD}_{y}^{2})\unicode[STIX]{x1D713}_{x}^{2}-2\unicode[STIX]{x1D6FD}_{x}\unicode[STIX]{x1D6FD}_{y}\unicode[STIX]{x1D713}_{x}\unicode[STIX]{x1D713}_{y}}{1+\unicode[STIX]{x1D6FD}_{x}^{2}+\unicode[STIX]{x1D6FD}_{y}^{2}}-\frac{\unicode[STIX]{x1D6FD}}{F^{2}}.\end{eqnarray}$$

Our numerical scheme presented in § 3 solves for the velocity potential $\unicode[STIX]{x1D713}$ . As such, the pressure field $p(x,y,\unicode[STIX]{x1D6FD})$ can be computed at very little additional cost and thus, in principle, the total force acting on the bottom topography can be calculated numerically by integrating the pressure over the bottom boundary.

To keep the scheme general, the bottom topography is left arbitrary in our formulation. However, for the purposes of presenting results, we make use of the bottom shape

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}(x,y)=-1+\unicode[STIX]{x1D716}\text{e}^{-((x-b_{1})^{2}+(y-b_{2})^{2})/2\unicode[STIX]{x1D6FF}^{2}}, & & \displaystyle\end{eqnarray}$$

which represents a Gaussian bump centred at $(b_{1},b_{2})$ with height $\unicode[STIX]{x1D716}$ and ‘width’ $\unicode[STIX]{x1D6FF}$ . We shall also use combinations of Gaussian bumps with different centres and widths.

2.2 Boundary integral formulation

In order to compute numerical solutions to (2.1) subject to (2.2)–(2.6), the problem is now reformulated into an integral equation using Green’s second theorem on the functions $\unicode[STIX]{x1D6F7}-x$ in the flow domain, with the Green’s function

(2.10) $$\begin{eqnarray}\frac{1}{R_{PQ}}=\frac{1}{4\unicode[STIX]{x03C0}}\frac{1}{((x-x^{\ast })^{2}+(y-y^{\ast })^{2}+(z-z^{\ast })^{2})^{1/2}}.\end{eqnarray}$$

Omitting the details, this reformulation results in the integral equation

(2.11) $$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x03C0}(\unicode[STIX]{x1D719}(x^{\ast },y^{\ast })-x^{\ast })=I_{1}+I_{2}+I_{3}+I_{4}, & & \displaystyle\end{eqnarray}$$

where

(2.12) $$\begin{eqnarray}\displaystyle I_{1} & = & \displaystyle \iint _{\mathbb{R}^{2}}(\unicode[STIX]{x1D719}(x,y)-\unicode[STIX]{x1D719}(x^{\ast },y^{\ast })-x+x^{\ast })\nonumber\\ \displaystyle & & \displaystyle \quad \times \,K_{1}(x,y,x^{\ast },y^{\ast },\unicode[STIX]{x1D701},\unicode[STIX]{x1D701}_{x},\unicode[STIX]{x1D701}_{y},\unicode[STIX]{x1D701}(x^{\ast },y^{\ast }))\,\text{d}x\,\text{d}y,\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle I_{2} & = & \displaystyle \iint _{\mathbb{R}^{2}}\unicode[STIX]{x1D701}_{x}(x,y)K_{2}(x,y,x^{\ast },y^{\ast },\unicode[STIX]{x1D701},\unicode[STIX]{x1D701}_{x},\unicode[STIX]{x1D701}_{y},\unicode[STIX]{x1D701}(x^{\ast },y^{\ast }))\,\text{d}x\,\text{d}y,\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle I_{3} & = & \displaystyle \iint _{\mathbb{R}^{2}}(\unicode[STIX]{x1D713}(x,y)-x)\nonumber\\ \displaystyle & & \displaystyle \times \,-\!K_{1}(x,y,x^{\ast },y^{\ast },\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FD}_{x},\unicode[STIX]{x1D6FD}_{y},\unicode[STIX]{x1D701}(x^{\ast },y^{\ast }))\,\text{d}x\,\text{d}y,\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle I_{4} & = & \displaystyle \iint _{\mathbb{R}^{2}}\unicode[STIX]{x1D6FD}_{x}(x,y)K_{2}(x,y,x^{\ast },y^{\ast },\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FD}_{x},\unicode[STIX]{x1D6FD}_{y},\unicode[STIX]{x1D701}(x^{\ast },y^{\ast }))\,\text{d}x\,\text{d}y,\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle K_{1}(x,y,x^{\ast },y^{\ast },a,b,c,d)=\frac{a-d-(x-x^{\ast })b-(y-y^{\ast })c}{((x-x^{\ast })^{2}+(y-y^{\ast })^{2}+(a-d)^{2})^{3/2}}, & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle K_{2}(x,y,x^{\ast },y^{\ast },a,b,c,d)=\frac{1}{\sqrt{(x-x^{\ast })^{2}+(y-y^{\ast })^{2}+(a-d)^{2}}}. & \displaystyle\end{eqnarray}$$

Very similar boundary integral formulations have been applied by other authors for steady three-dimensional flows – for example, for flows involving submerged singularities (Forbes Reference Forbes1989; Forbes & Hocking Reference Forbes and Hocking2005; Pethiyagoda et al. Reference Pethiyagoda, McCue, Moroney and Back2014a ,Reference Pethiyagoda, McCue and Moroney b ) and flows past pressure distributions (Părău & Vanden-Broeck Reference Părău and Vanden-Broeck2002, Reference Părău and Vanden-Broeck2011; Părău et al. Reference Părău, Vanden-Broeck and Cooker2005a ,Reference Părău, Vanden-Broeck and Cooker b , Reference Părău, Vanden-Broeck and Cooker2007a ,Reference Părău, Vanden-Broeck and Cooker b ). The key difference here is that we have allowed for an arbitrary bottom topography, which makes (2.11) more complicated.

In summary, we have reformulated our problem with the use of a boundary integral method so that, for a given bottom topography $\unicode[STIX]{x1D6FD}(x,y)$ and Froude number $F$ , we are required to solve the integral equation (2.11) and the dynamic condition (2.7) for the two velocity potentials $\unicode[STIX]{x1D719}(x,y)$ , $\unicode[STIX]{x1D713}(x,y)$ and the surface elevation $\unicode[STIX]{x1D701}(x,y)$ . For boundary conditions we impose the radiation condition (2.5)–(2.6).

3.1 Collocation scheme

The governing equations (2.7) and (2.11) are highly nonlinear. For the most gentle case in which the bottom topography is almost flat, we can simplify these equations to produce a linear problem, as detailed in § 4. However, to capture the full nonlinear features that arise for distinctively non-flat bottom topographies, we solve the equations numerically. With that in mind, the infinite plane is truncated such that $-\infty <x<\infty$ becomes $x_{1}<x<x_{N}$ and $-\infty <y<\infty$ becomes $y_{1}<y<y_{M}$ . Mesh points are then introduced as $x_{i}=x_{1}+(i-1)\unicode[STIX]{x0394}x$ , $i=1,2,\ldots ,N$ and $y_{j}=y_{1}+(j-1)\unicode[STIX]{x0394}y$ , $j=1,2,\ldots ,M$ . The point $x_{1}$ is chosen as a suitably far-upstream value such that the surface is considered flat there. The point $y_{1}$ is chosen such that $y_{1}=-y_{M}$ , although the scheme allows for asymmetrical meshes. The vector of $3(N+1)M$ unknowns $\boldsymbol{u}$  is

(3.1) $$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{u}= [\!\unicode[STIX]{x1D701}_{1,1},(\unicode[STIX]{x1D701}_{x})_{1,1},(\unicode[STIX]{x1D701}_{x})_{2,1},\ldots ,(\unicode[STIX]{x1D701}_{x})_{N,1},\unicode[STIX]{x1D701}_{1,2},(\unicode[STIX]{x1D701}_{x})_{1,2},\ldots ,(\unicode[STIX]{x1D701}_{x})_{N,M},\nonumber\\ \displaystyle & & \displaystyle \qquad \unicode[STIX]{x1D719}_{1,1},(\unicode[STIX]{x1D719}_{x})_{1,1},\ldots ,(\unicode[STIX]{x1D719}_{x})_{N,M},\unicode[STIX]{x1D713}_{1,1},(\unicode[STIX]{x1D713}_{x})_{1,1},\ldots ,(\unicode[STIX]{x1D713}_{x})_{N,M}\! ].\end{eqnarray}$$

The components of $\boldsymbol{u}$ are ordered by slicing the domain in the $x$ -direction. The values of the functions at the upstream truncation are placed before the derivatives of the corresponding slice.

The problem is solved using an inexact Newton method with an equation of the form $\boldsymbol{F}(\boldsymbol{u})=\mathbf{0}$ . First, the values for the equation for the free surface and the velocity potential on the free surface and the bottom surface are calculated using the trapezoidal rule. The singular integral $I_{2}$ can be written as $I_{2}=I_{2}^{\prime }+I_{2}^{\prime \prime }$ , where

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle I_{2}^{\prime }=\int _{y_{1}}^{y_{M}}\int _{x_{1}}^{x_{N}}(\unicode[STIX]{x1D701}_{x}(x,y)K_{2}-\unicode[STIX]{x1D701}_{x}(x^{\ast },y^{\ast })S_{2})\,\text{d}x\,\text{d}y, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle I_{2}^{\prime \prime }=\unicode[STIX]{x1D701}_{x}(x^{\ast },y^{\ast })\int _{y_{1}}^{y_{M}}\int _{x_{1}}^{x_{N}}S_{2}\,\text{d}x\,\text{d}y, & \displaystyle\end{eqnarray}$$

with

(3.4) $$\begin{eqnarray}S_{2}=\frac{1}{(A(x-x^{\ast })^{2}+B(x-x^{\ast })(y-y^{\ast })+C(y-y^{\ast })^{2})^{1/2}},\end{eqnarray}$$

and

(3.5a-c ) $$\begin{eqnarray}A=1+\unicode[STIX]{x1D701}_{x}^{2}(x^{\ast },y^{\ast }),\quad B=2\unicode[STIX]{x1D701}_{x}(x^{\ast },y^{\ast })\unicode[STIX]{x1D701}_{y}(x^{\ast },y^{\ast }),\quad C=1+\unicode[STIX]{x1D701}_{y}^{2}(x^{\ast },y^{\ast }).\end{eqnarray}$$

The integral $I_{2}^{\prime \prime }$ can be evaluated analytically to give

(3.6) $$\begin{eqnarray}\displaystyle \iint \frac{\text{d}s\,\text{d}t}{\sqrt{As^{2}+Bst+Ct^{2}}} & = & \displaystyle \frac{t}{\sqrt{A}}\ln (2As+Bt+2\sqrt{A(As^{2}+Bst+Ct^{2})})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{s}{\sqrt{C}}\ln (2Ct+Bs+2\sqrt{C(As^{2}+Bst+Ct^{2})}).\qquad\end{eqnarray}$$

Evaluating the boundary integral equation (2.11) at the half-mesh points $(x_{i+1/2},y_{j},i=1,2,\ldots ,N-1,j=1,2,\ldots ,M)$ produces $(N-1)M$ equations. Evaluating the boundary integral equation (2.11) at the half-mesh points on the bottom surface, that is, replacing $\unicode[STIX]{x1D701}(x^{\ast },y^{\ast })$ with $\unicode[STIX]{x1D6FD}(x^{\ast },y^{\ast })$ and $\unicode[STIX]{x1D719}(x^{\ast },y^{\ast })$ with $\unicode[STIX]{x1D713}(x^{\ast },y^{\ast })$ , produces another $(N-1)M$ equations. A further $(N-1)M$ equations come from evaluating the dynamic condition (2.7) at the half-mesh points. The last $6M$ equations come from the radiation condition outlined in Scullen (Reference Scullen1998),

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle x_{1}(\unicode[STIX]{x1D701}_{x})_{1,j}+n\unicode[STIX]{x1D701}_{1,j}=0, & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle x_{1}(\unicode[STIX]{x1D701}_{xx})_{1,j}+n(\unicode[STIX]{x1D701}_{x})_{1,j}=0, & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle x_{1}((\unicode[STIX]{x1D719}_{x})_{1,j}-1)+n(\unicode[STIX]{x1D719}_{1,j}-x)=0, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle x_{1}(\unicode[STIX]{x1D719}_{xx})_{1,j}+n((\unicode[STIX]{x1D719}_{x})_{1,j}-1)=0, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle x_{1}((\unicode[STIX]{x1D713}_{x})_{1,j}-1)+n(\unicode[STIX]{x1D713}_{1,j}-x)=0, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle x_{1}(\unicode[STIX]{x1D713}_{xx})_{1,j}+n((\unicode[STIX]{x1D713}_{x})_{1,j}-1)=0, & \displaystyle\end{eqnarray}$$

for $j=1,2,\ldots ,M$ . The value $n=0.05$ represents the upstream decay rate; this particular value was also used by Pethiyagoda et al. (Reference Pethiyagoda, McCue, Moroney and Back2014a ). We now have $3(N+1)M$ nonlinear equations (that make up the vector $\boldsymbol{F}$ ) for the $3(N+1)M$ unknowns in (3.1).

3.2 Jacobian-free Newton–Krylov method

A Jacobian-free Newton–Krylov (JFNK) method is used to solve the system

(3.13) $$\begin{eqnarray}\boldsymbol{F}(\boldsymbol{u})=\mathbf{0}.\end{eqnarray}$$

This method uses a damped Newton iteration

(3.14) $$\begin{eqnarray}\boldsymbol{u}_{k+1}=\boldsymbol{u}_{k}+\unicode[STIX]{x1D707}_{k}\unicode[STIX]{x1D739}\boldsymbol{u}_{k},\end{eqnarray}$$

with a simple line search to find $\unicode[STIX]{x1D707}_{k}$ such that there is a sufficient decrease in the nonlinear residual. The Newton step is found by satisfying

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D645}(\boldsymbol{u}_{k})\unicode[STIX]{x1D739}\boldsymbol{u}_{k}=-\boldsymbol{F}(\boldsymbol{u}_{k}),\end{eqnarray}$$

where $\unicode[STIX]{x1D645}$ is the Jacobian matrix. A generalised minimum residual (GMRES) method (Saad & Schultz Reference Saad and Schultz1986) with right preconditioning is used to solve for $\unicode[STIX]{x1D739}\boldsymbol{u}_{k}$ . The approximate solution for $\unicode[STIX]{x1D739}\boldsymbol{u}_{k}$ is found, after $m$ iterations, by projecting onto the Krylov subspace

(3.16) $$\begin{eqnarray}{\mathcal{K}}_{m}(\unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1},\boldsymbol{F}_{k})=\text{span}\{\boldsymbol{F}_{k},\unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1}\boldsymbol{F}_{k},\ldots ,(\unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1})^{m-1}\boldsymbol{F}_{k}\},\end{eqnarray}$$

where $\unicode[STIX]{x1D645}_{k}=\unicode[STIX]{x1D645}(\boldsymbol{u}_{k})$ and $\boldsymbol{F}_{k}=\boldsymbol{F}(\boldsymbol{u}_{k})$ . The preconditioner matrix $\unicode[STIX]{x1D64B}\approx \unicode[STIX]{x1D645}_{k}$ is a sparse approximation to the full Jacobian $\unicode[STIX]{x1D645}_{k}$ and will be discussed further later. The purpose of the preconditioner is to reduce the number of iterations of the GMRES algorithm by reducing the dimension of the Krylov subspace necessary to find an accurate value for $\unicode[STIX]{x1D739}\boldsymbol{u}_{k}$ .

The benefit of using a Krylov subspace method is that it does not require the full Jacobian to be formed. The Jacobian is only actioned on vectors to form the basis of the preconditioned Krylov subspace ${\mathcal{K}}_{m}$ . First-order quotients can be used to approximate these actions without having to explicitly form $\unicode[STIX]{x1D645}_{k}$ :

(3.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1}\boldsymbol{v}\approx \frac{\boldsymbol{F}(\boldsymbol{u}_{k}+h\unicode[STIX]{x1D64B}^{-1}\boldsymbol{v})-\boldsymbol{F}(\boldsymbol{u}_{k})}{h}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{v}$ is an arbitrary vector used to form the Krylov subspace and $h$ is a sufficiently small shift (Brown & Saad Reference Brown and Saad1990). Due to both the action of the Jacobian and the Newton correction both being approximated, we are left with an inexact Newton method, which leads to only superlinear and not quadratic convergence (Knoll & Keyes Reference Knoll and Keyes2004).

The aim when forming the preconditioner matrix $\unicode[STIX]{x1D64B}$ is to construct an approximation to the Jacobian $\unicode[STIX]{x1D645}_{k}$ that is easy to form and factorise, and exhibits a clustering of eigenvalues in the spectrum of the preconditioned Jacobian $\unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1}$ . A good starting point for such a preconditioner is to consider a simplified form of the equations. In our case, we will be applying our integral formulation to the linearised problem that arises in the limit as the bottom topography approaches a flat bottom.

4.1 Exact solution

In this section we consider the linearised version of (2.7) and (2.11) which arises from taking the limit that the perturbation of the bottom topography from a flat surface vanishes. This is done mathematically via a perturbation expansion in a parameter $\unicode[STIX]{x1D716}$ which is a measure of the ‘height’ of the bottom disturbance. For the particular example (2.9), the parameter $\unicode[STIX]{x1D716}$ is defined as the height of the Gaussian.

By taking the linear limit, the conditions on $z=\unicode[STIX]{x1D701}(x,y)$ and $z=\unicode[STIX]{x1D6FD}(x,y)$ are projected onto the planes $z=0$ and $z=-1$ , respectively. Omitting the details, the linearised problem is to solve Laplace’s equation

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6F7}=0,\quad -1<z<0,\end{eqnarray}$$

subject to the kinematic, dynamic and radiation conditions (2.2)–(2.5), which are linearised under the perturbation expansion to become

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}_{x}=\unicode[STIX]{x1D6F7}_{z}\quad z=0, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{x}=\unicode[STIX]{x1D6F7}_{z}\quad z=-1, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}_{x}-1+\frac{\unicode[STIX]{x1D701}}{F^{2}}=0\quad z=0, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6F7}_{x},\unicode[STIX]{x1D6F7}_{y},\unicode[STIX]{x1D6F7}_{z})\rightarrow (1,0,0)\quad x\rightarrow -\infty , & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}\rightarrow 0\quad x\rightarrow -\infty . & \displaystyle\end{eqnarray}$$

The linear problem (4.1)–(4.6) can be solved exactly using Fourier transforms to give

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}(x,y) & = & \displaystyle \frac{F^{2}}{2\unicode[STIX]{x03C0}^{2}}\int _{-\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x03C0}/2}\int _{0}^{\infty }\frac{k^{2}\,\text{sech}(k)\tilde{\unicode[STIX]{x1D6FD}}(k,\unicode[STIX]{x1D713})\cos (k[|x|\cos \unicode[STIX]{x1D713}+y\sin \unicode[STIX]{x1D713}])}{kF^{2}-\sec ^{2}\unicode[STIX]{x1D713}\tanh k}\,\text{d}k\,\text{d}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{2F^{2}H(x)}{\unicode[STIX]{x03C0}}\int _{\unicode[STIX]{x1D713}_{0}}^{\unicode[STIX]{x03C0}/2}\frac{k_{1}^{2}\,\text{sech}(k_{1})\tilde{\unicode[STIX]{x1D6FD}}(k_{1},\unicode[STIX]{x1D713})\sin (k_{1}x\cos \unicode[STIX]{x1D713})\cos (k_{1}y\sin \unicode[STIX]{x1D713})}{F^{2}-\sec ^{2}\unicode[STIX]{x1D713}\,\text{sech}^{2}k_{1}}\,\text{d}\unicode[STIX]{x1D713},\qquad\end{eqnarray}$$

where $H(x)$ is the Heaviside function, $\tilde{\unicode[STIX]{x1D6FD}}(k,\unicode[STIX]{x1D713})$ is the Fourier transform of $\unicode[STIX]{x1D6FD}(x,y)+1$ in polar coordinates, $\unicode[STIX]{x1D713}_{0}=0$ for $F<1$ , $\unicode[STIX]{x1D713}_{0}=\arccos (1/F)$ for $F>1$ and the path of $k$ -integration is taken below the pole $k=k_{1}$ , where $k_{1}$ is the real positive root of $kF^{2}-\sec ^{2}\unicode[STIX]{x1D713}\tanh k=0$ . This formula (4.7) provides an extremely good approximation to the free surface for flow regimes in which the bottom topography is almost flat.

It is interesting to note that the exact solution (4.7) looks very similar to that for free-surface flow past a pressure distribution $p(x,y)$ that has been applied to the surface (Pethiyagoda, McCue & Moroney Reference Pethiyagoda, McCue and Moroney2015). In fact, the two problems are closely related. Suppose we have a linear solution for flow past a bottom obstruction for a given bottom shape $z=\unicode[STIX]{x1D6FD}(x,y)$ and corresponding Fourier transform $\tilde{\unicode[STIX]{x1D6FD}}(k,\unicode[STIX]{x1D713})$ . Then for the related linear problem of flow past a pressure distribution, if the pressure $p(x,y)$ is chosen so that its Fourier transform is given by

(4.8) $$\begin{eqnarray}\tilde{p}(k,\unicode[STIX]{x1D713})=\frac{\text{sech}(k)\tilde{\unicode[STIX]{x1D6FD}}(k,\unicode[STIX]{x1D713})}{\unicode[STIX]{x1D716}F^{2}},\end{eqnarray}$$

then the wave pattern behind the pressure distribution is given by (4.7), except there is an addition term $-p(x,y)$ on the right-hand side. That is, in the linear regime, the wave pattern downstream from the disturbance is the same. Therefore, given that the problem of flow past an applied pressure is used as a model for studying ship wakes, there is a direct relationship between the wave patterns considered in this paper and ship waves. We comment further on this relationship in § 6.

4.2 Boundary integral formulation

Apart from benefiting from writing out the exact solution (4.7), the other reason for us to pursue the linearised problem is to construct a preconditioner for the fully nonlinear problem. For this purpose we reformulate the linear problem (4.1)–(4.6) using Green’s second theorem to give

(4.9) $$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x03C0}(\unicode[STIX]{x1D719}^{\ast }-x^{\ast }) & = & \displaystyle \iint _{\mathbb{R}^{2}}\frac{\unicode[STIX]{x1D701}_{x}}{\sqrt{(x-x^{\ast })^{2}+(y-y^{\ast })^{2}}}\,\text{d}x\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle +\,\iint _{\mathbb{R}^{2}}\frac{\unicode[STIX]{x1D6FD}_{x}}{\sqrt{(x-x^{\ast })^{2}+(y-y^{\ast })^{2}+1}}\,\text{d}x\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle +\,\iint _{\mathbb{R}^{2}}\frac{(\unicode[STIX]{x1D713}-x^{\ast })}{((x-x^{\ast })^{2}+(y-y^{\ast })^{2}+1)^{3/2}}\,\text{d}x\,\text{d}y,\end{eqnarray}$$

which is coupled to the dynamic condition

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{x}^{\ast }+\frac{\unicode[STIX]{x1D701}^{\ast }}{F^{2}}-1=0.\end{eqnarray}$$

By applying the same discretisation to that described in § 3.1, we can derive a completely analogous system of $3(N+1)M$ equations for the same $3(N+1)M$ unknowns (3.1). The resulting equations in their exact form are detailed in the supplementary material available at https://doi.org/10.1017/jfm.2018.254. We are interested in exploiting the structure and entries of the Jacobian, $\tilde{\unicode[STIX]{x1D645}}$ , of this new linear system.

Figure 2 shows the similarity between the nonlinear Jacobian $\unicode[STIX]{x1D645}_{k}$ and the Jacobian of the corresponding linear problem $\tilde{\unicode[STIX]{x1D645}}$ . The visualisation shows that $\tilde{\unicode[STIX]{x1D645}}$ shares the structure of $\unicode[STIX]{x1D645}_{k}$ , while being much easier to form as analytical formulae exist for the derivatives which are provided in their exact form in the supplementary material. Since analytical formulae exist for the elements of $\tilde{\unicode[STIX]{x1D645}}$ , the elements can be calculated independently and utilise parallel processing. We are also able to exploit the block structure of the matrix by forming the blocks separately. This is useful as $\tilde{\unicode[STIX]{x1D645}}_{2,3}=\tilde{\unicode[STIX]{x1D645}}_{3,2}$ , where $\tilde{\unicode[STIX]{x1D645}}$ is treated as a $3\times 3$ block matrix, and $\tilde{\unicode[STIX]{x1D645}}_{2,2}=\tilde{\unicode[STIX]{x1D645}}_{3,3}$ , which allows us to save time by not explicitly forming two blocks.

Figure 2.

Visualisation of the nonlinear and linear Jacobians showing the logarithm (base 10) of the magnitudes of the entries. Entries with larger magnitudes are represented by darker shades and zero entries are coloured white. Jacobians were generated from the problem on a $31\times 31$ mesh with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0.8$ , $F=0.6$ , and a single Gaussian bump on the bottom surface with scale parameter $\unicode[STIX]{x1D716}=0.1$ . (a) The full Jacobian for the completely nonlinear problem, which is generated using finite differences. (b) The Jacobian of the linearised problem. The linear Jacobian can be generated analytically by taking derivatives of the linearised equations. The two Jacobians clearly share a similar structure, with dense blocks and sparse blocks in the same places. The key differences are the top-middle block having more of a diagonal structure and the bottom-right block only having the block diagonal structure in the linear case. The structural similarities suggest that the linear Jacobian may be a good approximation for the nonlinear Jacobian.

An effective preconditioner has to reduce the dimension of the Krylov subspace required to find a solution. To test whether the linear Jacobian is an effective preconditioner, we can compare the eigenvalues, $\unicode[STIX]{x1D706}$ , of the nonlinear Jacobian and the preconditioned nonlinear Jacobian. Figure 3(a) shows the eigenvalues of the nonlinear Jacobian $\unicode[STIX]{x1D645}_{k}$ for a flat surface. Figure 3(b) shows the eigenvalues for the nonlinear Jacobian preconditioned with the linear Jacobian – that is, $\unicode[STIX]{x1D645}_{k}\tilde{\unicode[STIX]{x1D645}}^{-1}$ . The clustering of eigenvalues is much tighter than the non-preconditioned Jacobian. Thus, we chose that $\unicode[STIX]{x1D64B}=\tilde{\unicode[STIX]{x1D645}}$ .

Figure 3.

Eigenvalue plots of (a) the full nonlinear Jacobian $\unicode[STIX]{x1D645}_{k}$ for an initial guess of a flat surface with $F=0.6$ and $\unicode[STIX]{x1D716}=0.1$ on a $31\times 31$ mesh with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0.8$ and the full nonlinear Jacobian preconditioned with (b) the linear preconditioner $\unicode[STIX]{x1D645}_{k}\tilde{\unicode[STIX]{x1D645}}^{-1}$ and (c) the banded linear preconditioner $\unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1}$ . The effect of preconditioning is a tight clustering of the eigenvalues, $\unicode[STIX]{x1D706}$ , around $1$ . This tight clustering is desirable as a smaller spectrum of the preconditioned Jacobian, $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D645}_{k}\tilde{\unicode[STIX]{x1D645}}^{-1})$ , will result in a smaller Krylov subspace required to find a solution to the Newton step $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{k}$ .

Now that we have a preconditioner $\unicode[STIX]{x1D64B}$ that is cheap to form, we want to find a way to reduce the large memory requirement to store $\unicode[STIX]{x1D64B}$ . In the four dense blocks in $\unicode[STIX]{x1D64B}$ the magnitude of the elements decays to zero away from the main diagonal. This structure suggests that it is possible to store only the elements near the diagonal using a banded approximation of the dense blocks. Indeed, we follow this approach and use a banded version of $\tilde{\unicode[STIX]{x1D645}}$ for our preconditioner $\unicode[STIX]{x1D64B}$ . The details are described in the Appendix. Figure 3(c) shows the eigenvalues for the nonlinear Jacobian preconditioned with the banded preconditioner, $\unicode[STIX]{x1D645}_{k}\unicode[STIX]{x1D64B}^{-1}$ . While the clustering is not as tight as when using the full linear Jacobian, the dimension of the Krylov subspace required to find a solution is still significantly reduced.

5.1 Subcritical flow

Subcritical flows are defined by $F<1$ . In this regime, the wave pattern is characterised by the presence of both divergent and transverse waves. In figure 4, we provide some representative results for a fixed subcritical Froude number $F=0.6$ . In all cases, the bottom topography takes the form (2.9) with $\unicode[STIX]{x1D6FF}=0.5$ , which represents a single Gaussian-type bump with height $\unicode[STIX]{x1D716}$ on an otherwise flat bottom.

Figure 4(a,b) shows the wave pattern for the exact linear solution (4.7), which is valid for small bump heights, $\unicode[STIX]{x1D716}\ll 1$ . We can see that, for this linear solution, the wave pattern appears to be dominated by the transverse waves which run roughly perpendicular to the direction of flow. Similar free-surface profiles are provided in figure 4(c,d). This time the bump height was fixed to be $\unicode[STIX]{x1D716}=0.1$ and the solution was computed using our fully nonlinear numerical scheme. Thus we can see for small ( $\unicode[STIX]{x1D716}\ll 1$ ) to moderately small ( $\unicode[STIX]{x1D716}=0.1$ ) obstructions, the linear and nonlinear results are not significantly different.

For larger values of $\unicode[STIX]{x1D716}$ , differences emerge between the linear and nonlinear regimes, and these differences become more pronounced as $\unicode[STIX]{x1D716}$ increases. For example, in figure 4(e,f) we present results for $\unicode[STIX]{x1D716}=0.285$ , which represents a highly nonlinear solution. In figure 4(e), in particular, we see a remarkably nonlinear wave pattern with sharp divergent waves shaped like dorsal fins, whose amplitude is much higher than the transverse waves (these are reminiscent of the highly nonlinear solutions computed by Pethiyagoda et al. (Reference Pethiyagoda, McCue and Moroney2014b ) for flow past a submerged source singularity). It is important to emphasise that the distinctive shape of the sharp divergent waves are not predicted by linear theory, demonstrating that for sufficiently large disturbances to the free stream, a nonlinear formulation is required to accurately describe the wave pattern.

As an indication of the relative vertical scales involved, we present in figure 4(g) a plot of the two nonlinear surfaces along the centreline $y=0$ . Clearly the amplitude of the waves for $\unicode[STIX]{x1D716}=0.285$ is significantly higher than that for $\unicode[STIX]{x1D716}=0.1$ , as expected. Further, the waves for $\unicode[STIX]{x1D716}=0.285$ appear to have broader troughs and sharper crests. Another feature of figure 4(g) is that even though the Froude number is fixed, the wavelength is smaller for the more nonlinear solution. These are all commonly observed properties of nonlinear free-surface flows in two dimensions (Forbes & Schwartz Reference Forbes and Schwartz1982; Forbes Reference Forbes1985; McCue & Forbes Reference McCue and Forbes2002).

We note that for $F=0.6$ , we did not compute solutions for $\unicode[STIX]{x1D716}>0.285$ . From a mathematical perspective, we expect that solutions will exist up to some maximum value of $\unicode[STIX]{x1D716}$ , at which point the surface will be in some limiting configuration which is likely to have the maximum surface height at $F^{2}/2$ (Bernoulli’s equation (2.4) does not allow for the free surface to be higher than $F^{2}/2$ ). This limiting configuration would be some type of complicated three-dimensional analogue of Stokes’ limiting configuration in two dimensions (with a $120^{\circ }$ corner at the crest).

Returning to the free-surface profile in figure 4(c,d), small-amplitude waves appear ahead of the bottom disturbance which are not physical and are caused by truncating the domain upstream. These spurious numerical artefacts appear in all of our solutions for subcritical Froude numbers, and can be moderated by truncating further upstream (which, for a fixed grid spacing, involves increasing the number of grid points) or by varying the parameter $n$ in the radiation condition (3.7)–(3.12). Note that these numerical waves appear larger in figure 4(c,d) than in figure 4(e,f); however, in fact, they are similar in size (see figure 4 g, for example). It is the vertical scaling of the surfaces that causes this illusion.

Figure 5.

Numerical solutions for supercritical flow past a single bump (2.9) centred at the origin. Solutions were computed for $-21\leqslant x\leqslant 51$ , $-36\leqslant y\leqslant 36$ using a $281\times 281$ mesh with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0.257$ , $x_{1}=-21$ , $\unicode[STIX]{x1D6FF}=3$ and $F=3$ . (a,b) The linear solution. (c,d) The nonlinear solution with $\unicode[STIX]{x1D716}=0.1$ . (e,f) The nonlinear solution with $\unicode[STIX]{x1D716}=2.75$ . As an indication of the relative vertical scales in these plots, the maximum height for the solution in (c,d) is $\unicode[STIX]{x1D701}(x,0)=1.091$ while the maximum height for the solution in (e,f) is $\unicode[STIX]{x1D701}(x,0)=2.556$ . (g) A cut along the centreline for the solution with $\unicode[STIX]{x1D716}=2.75$ , including the free surface $\unicode[STIX]{x1D701}(x,0)$ and the bottom $\unicode[STIX]{x1D6FD}(x,0)$ .

5.2 Supercritical flow

Supercritical flows are for $F>1$ . In figure 5, we present supercritical solutions for the representative value $F=3$ . Again, we have chosen the bottom obstruction to be a single bump of the form (2.9), but now with $\unicode[STIX]{x1D6FF}=3$ . Figure 5(a,b) shows free-surface profiles for the exact linear solution (4.7), while figure 5(c,d) shows the surfaces for the bump height $\unicode[STIX]{x1D716}=0.1$ . The first point to note here is that the wave pattern appears very different from that presented for subcritical flows. This qualitative difference is well known, and is for the most part due to supercritical wave patterns being entirely made up of divergent waves. According to linear theory, the mathematical reason for the absence of transverse waves is that the dispersion relation $kF^{2}-\sec ^{2}\unicode[STIX]{x1D713}\tanh k=0$ has no real positive roots for $0\leqslant \unicode[STIX]{x1D713}<\arccos (1/F)$ . Physically, if we consider an extremely localised disturbance, then for supercritical flow the surface elevation at a single point downstream is due to a wave that began propagating out from the disturbance at the group velocity some time ago, while for subcritical flow the elevation at a single point is due to two waves that began propagating at two different previous times (geometric constructions like these are explained in detail in Lamb (Reference Lamb1932), Wehausen & Laitone (Reference Wehausen and Laitone1960), Soomere (Reference Soomere2007), Pethiyagoda et al. (Reference Pethiyagoda, Moroney, MacFarlane, Binns and McCue2018a ), for example). The second point to note is that for such a small value as $\unicode[STIX]{x1D716}=0.1$ , the nonlinear wave pattern appears very similar to the linear solution.

We are motivated to increase the bump height $\unicode[STIX]{x1D716}$ to observe the effects of nonlinearity. It turns out that we are able to compute numerical solutions to the fully nonlinear problem for much larger values of $\unicode[STIX]{x1D716}$ . For example, in figure 5(e,f) we show free-surface profiles for $\unicode[STIX]{x1D716}=2.75$ . Given that $\unicode[STIX]{x1D716}$ measures the height of the Gaussian (2.9), this value of $\unicode[STIX]{x1D716}$ corresponds to a significant deviation from a flat bottom. To illustrate the scales involved, we have presented in figure 5(g) a slice of the solution along the centreline $y=0$ , including the free surface $z=\unicode[STIX]{x1D701}(x,0)$ and the bottom boundary $z=\unicode[STIX]{x1D6FD}(x,0)$ . This type of solution is very interesting because the bottom disturbance is much higher than the depth of the channel far upstream, even though the bump itself does not pierce the surface.

Further, we can see some features in figure 5(e,f) that are not present in the linear solution. For example, the V-shaped wave ridge is much steeper for $\unicode[STIX]{x1D716}=2.75$ than the linear solution; the angle of this V-shaped ridge (the apparent wake angle) is larger than that for the linear solution, and there appears to be a second, smaller, V-shaped wedge that forms a distance behind the disturbance. This is a new feature that is not able to be predicted by linear theory.

To investigate these nonlinear effects further, we show in figure 6(a) the dependence of the apparent wake angle on the obstruction height $\unicode[STIX]{x1D716}$ for two supercritical Froude numbers, $F=3$ and $4$ . These two plots show how the wake angle is significantly lower for the larger Froude number, which is true for both linear and nonlinear regimes. In addition, the increase in wake angle with nonlinearity is clear from these plots. Also shown in figure 6(b) is a plot of the maximum steepness of the surface versus the obstruction height $\unicode[STIX]{x1D716}$ . Here wave steepness is defined to be the maximum slope on the outside of the (main) V-shaped ridge that extends to the far field. We see that the wave steepness clearly increases with nonlinearity. Presumably this steepness will continue to increase as $\unicode[STIX]{x1D716}$ increases, until the surface near the origin reaches its limiting height of $F^{2}/2$ .

Figure 6.

Plots of (a) the apparent wake angle and (b) the wave steepness against bump height $\unicode[STIX]{x1D716}$ for $F=3$ (solid blue circles) and $4$ (open red circles). The apparent wake angle is given by the angle between the V-shaped ridge and the centreline. The wave steepness is defined as the maximum slope of the surface on the outside of the V-shaped ridge. These results demonstrate that nonlinearity has a significant effect on the wave pattern for supercritical flow regimes.

5.3 Multiple bumps

An advantage of formulating our numerical scheme with an arbitrarily defined bottom topography $z=\unicode[STIX]{x1D6FD}(x,y)$ is that solutions can easily be computed for flow over multiple bumps on the bottom surface. Two such flow configurations are shown in figure 7 for a fixed Froude number $F=0.6$ . Figure 7(a,b) shows the numerical solution to the nonlinear problem with two bumps positioned symmetrically across the $x$ -axis. The wave trains appear independent for sufficiently small $x$ and then interfere with each other at approximately $x=30$ . In figure 7(c,d) we show the wave pattern for flow past two bumps positioned asymmetrically across the $x$ -axis. In this configuration, one of the bumps is further downstream than the other one. Again, the wave trains interfere with each other, as expected. Note for this configuration there is no direct analogue in two dimensions.

We now consider the situation in which there are two bumps positioned inline along the $x$ -axis. Four free-surface profiles for this configuration are presented in figure 8, again for $F=0.6$ . The profile in figure 8(a,b) is a linear solution computed for two bumps that are positioned 4 wavelengths apart, giving rise to constructive interference of the transverse waves. In contrast, the profile in figure 8(c,d) is a linear solution computed for two bumps that are positioned 4.5 wavelengths apart, leading to destructive interference. Here a single wavelength is given by $2\unicode[STIX]{x03C0}/k$ , where $k$ is the real positive root of the linear dispersion relation $kF^{2}-\tanh k=0$ (for $F=0.6$ this linear wavelength is approximately $2.28$ ). As the transverse waves decay in amplitude (like $x^{-1/2}$ for $x\gg 1$ ) the interference is not absolute (for example, the waves do not completely cancel each other out in the new figure 8 c,d). However, the effects are strong and, for the case of constructive interference, we see the resulting transverse waves downstream become larger than the divergent waves, which is unusual for this Froude number.

Analogous nonlinear profiles are presented in figure 8(e,f) and figure 8(g,h). In this nonlinear case there is no obvious wavelength as the linear theory no longer applies. Instead, we have computed these solutions for bumps that are approximately 4 and 4.5 wavelengths apart, where a single wavelength is estimated from a nonlinear solution with a single hump. Again, the visual differences between the two cases are clear from these figures. To support these results, in figure 9 we have presented centreline plots which more clearly show the effects of (transverse) wave interference. Similar plots can be generated for bumps that are separated by any number of wavelengths, and the wave profiles exhibit similar properties. These results are reminiscent of the wave interference effects that occur at the stern of a steadily moving ship, although in that context it is also of interest to study interference by divergent waves (Noblesse et al. Reference Noblesse, He, Zhu, Hong, Zhang, Zhu and Yang2014; Zhang et al. Reference Zhang, He, Zhu, Yang, Li, Zhu, Lin and Noblesse2015; Zhu et al. Reference Zhu, He, Zhang, Wu, Wan, Zhu and Noblesse2015).

Figure 7.

Numerical solutions for subcritical flow past two bumps of the form (2.9). Solutions were computed for $-7\leqslant x\leqslant 53$ , $-30\leqslant y\leqslant 30$ using a $281\times 281$ mesh with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0.2143$ , $F=0.6$ , bump height $\unicode[STIX]{x1D716}=0.2$ , $\unicode[STIX]{x1D6FF}=0.5$ and the upstream truncation point $x_{1}=-7$ . (a,b) The solution for flow past two bumps positioned symmetrically across the $x$ -axis at $(0,8)$ and $(0,-8)$ . (c,d) The solution for flow past two bumps offset in the $x$ -direction at $(10,10)$ and $(0,-10)$ .

Figure 8.

Free-surface profiles for linear and nonlinear flow past two submerged bumps positioned inline along the $x$ -axis, calculated for $\unicode[STIX]{x1D716}=0.1$ , $\unicode[STIX]{x1D6FF}=0.5$ and $F=0.6$ . The separation distance between the humps is chosen so that the wave profiles exhibit constructive interference (linear in a,b and nonlinear in e,f) and destructive interference (linear in c,d and nonlinear in g,h) of the transverse waves. The nonlinear solutions were generated for $-7\leqslant x\leqslant 53$ using a $281\times 281$ mesh with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0.2143$ . Note the centreline plots for these solutions are provided in Figure 9.

Figure 9.

Centreline plots exhibiting both destructive (blue solid) and constructive (red dashed) interference along the $y$ -axis, computed for $\unicode[STIX]{x1D716}=0.1$ , $\unicode[STIX]{x1D6FF}=0.5$ and $F=0.6$ , as in Figure 8. (a) Linear solutions. (b) Nonlinear solutions.

5.4 Flow past a crater

We now turn our attention to flow past an isolated crater, which is easy to implement with our numerical scheme (see Zhang & Zhu (Reference Zhang and Zhu1996a ), Wade et al. (Reference Wade, Binder, Mattner and Denier2017) for the two-dimensional analogue). In particular, we consider flow over an inverted Gaussian, which is (2.9) with $\unicode[STIX]{x1D716}<0$ . Figure 10(a,b) shows surface profiles for $\unicode[STIX]{x1D716}=-0.1$ , which is a relatively shallow crater, while figure 10(c,d) is for $\unicode[STIX]{x1D716}=-0.62$ , a relatively deep crater. The solution for $\unicode[STIX]{x1D716}=-0.1$ has a wave pattern that is very similar to the linear solution (not shown) and is also an almost-exact reflection of the surface for $\unicode[STIX]{x1D716}=0.1$ (see figure 4 c,d). On the other hand, the highly nonlinear solution for $\unicode[STIX]{x1D716}=-0.62$ is interesting because its wave pattern appears very different from both the linear solution and the nonlinear solution for flow over a standard Gaussian. For example, even though this is a highly nonlinear solution, the divergent waves do not resemble the steep dorsal fins we see in figure 4(e).

Figure 10.

Nonlinear solutions for subcritical flow past a submerged crater centred at the origin. Solutions were computed using a $281\times 281$ mesh with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0.0857$ , $x_{1}=-7$ , $\unicode[STIX]{x1D6FF}=0.5$ , $F=0.6$ with $\unicode[STIX]{x1D716}=-0.1$ in (a,b) and $\unicode[STIX]{x1D716}=-0.62$ in (c,d). The solutions show that for the near-linear case (a,b), the wave patterns are the inverse of the solution for flow past a bump; however, for the more nonlinear case (c,d), the first wave crest is much higher than the rest of the waves and the transverse waves appear to be out of phase with the divergent waves. As an indication of the relative vertical scales in these plots, the maximum height for the solution in (a,b) is 0.020 while the maximum height for the solution in (c,d) is 0.073.