2.1 Tulin’s water-wave model

The essential derivation behind Tulin’s model begins from Bernoulli’s equation applied to a free surface with streamline, $\unicode[STIX]{x1D713}=0$ ,

where $q$ is the fluid speed, $\unicode[STIX]{x1D703}$ the streamline angle, $\unicode[STIX]{x1D719}$ the potential, $\unicode[STIX]{x1D716}$ is given by (1.2). If the sinusoidal term is split according to the identity

then (2.1) can be written in complex-valued form as

where $G=(q\text{e}^{-\text{i}\unicode[STIX]{x1D703}})^{3}$ is an analytic function of the complex potential, $w=\unicode[STIX]{x1D719}+\text{i}\unicode[STIX]{x1D713}$ , and the above is evaluated on $\unicode[STIX]{x1D713}=0$ where $\text{Re}({\mathcal{P}})=4\sin ^{3}\unicode[STIX]{x1D703}/q^{3}$ .

The rather curious substitution of (2.2) is attributed to Davies (Reference Davies1951, p. 478), who had argued that if ${\mathcal{P}}$ is considered to be small, then (2.3) yields a linearized version of Bernoulli’s equation (in $G$ ) that preserves the essential nonlinearity (in $q^{3}$ ), which describes the structure of a steep gravity wave with a crest angle near $120^{\circ }$ . As explained by Davies, the rationale is that the nonlinear variation in the speed of a sharp-crested wave is an essential component to preserve, but not so much the variation in $\unicode[STIX]{x1D703}$ .

Inspired by this idea, Tulin considered the extension to general free-surface flows over a moving body. Since the function in the curly braces of (2.3) is an analytic function in $w$ , except at isolated points in the complex plane, then analytic continuation implies

where by matching to uniform flow, $G\rightarrow 1$ and thus ${\mathcal{Q}}\rightarrow \text{i}$ in the upstream limit. The function ${\mathcal{Q}}(w)$ is purely imaginary on $\unicode[STIX]{x1D713}=0$ . In the case of flow without a body or bottom boundary, ${\mathcal{Q}}\equiv \text{i}$ , but otherwise, ${\mathcal{Q}}$ will encode the effect of the obstructions through its singular behaviour in the complex plane. If the nonlinear contribution, ${\mathcal{P}}$ , is neglected in (1.1), then the solution can be written as

where the lower limit of integration should be chosen to zero the constant of integration. Hence (2.5) yields Tulin’s exact theory.

This process seems to yield the exponentially small surface waves at low speeds. However, Tulin noted that as $\unicode[STIX]{x1D716}\rightarrow 0$ , there would be locations on the free surface where ${\mathcal{Q}}=0$ , and this would lead to unbounded steepness. He wrote of

He thus proposed the following result:

The matter was left at that, and in the three decades following Tulin’s ingenious paper, the peculiarities surrounding the asymptotic breakdown of (1.1) were never directly readdressed (although we mention the paper by Vanden-Broeck & Miloh (Reference Vanden-Broeck and Miloh1995) which develops numerical solutions for the case ${\mathcal{Q}}=\text{const.}$ and ${\mathcal{P}}\neq 0$ ).

Note that in addition to the investigation in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ , Tulin had also intended to produce an exact reduction, where the ${\mathcal{P}}$ term in (1.1) was handled through a theoretically posited nonlinear coordinate transformation. However, it is never clear how this transformation is used or derived in practice, except through nonlinear computation of the full solutions.

2.2 Tuck’s linearized integral equation

Independently from Tulin, Tuck later presented a series of papers (Tuck Reference Tuck1990, Reference Tuck1991a ,Reference Tuck b ) where he attempted to distil the wave-making properties of wave–body problems into a linear singular equation. The equation presented (equation (22) of Tuck (Reference Tuck1991a )) was

where $y(\unicode[STIX]{x1D719})$ is the height of the free surface, ${\mathcal{F}}(\unicode[STIX]{x1D719})$ is a function related to the moving body, and $\mathscr{H}$ is an integral operator known as the Hilbert transform (to be introduced in § 3). The difficulty in solving (2.6) is that the Hilbert transform is a global operator, requiring values of $y$ over the entire domain.

In the Reference Tuck1991a work, Tuck explained how the action of $\mathscr{H}$ could be viewed as similar to that of the differential operator $\text{d}/\text{d}\unicode[STIX]{x1D719}$ . This reduction was mostly pedagogic in nature, but Tuck was motivated by the fact that $\mathscr{H}$ and the differential operator will act similarly in the case of sinusoidal functions. By replacing $\mathscr{H}$ by $\unicode[STIX]{x1D716}\text{d}/\text{d}\unicode[STIX]{x1D719}$ , he went on to study the various properties of the singular differential equation

that depend on the specification of the ‘body’ given by ${\mathcal{F}}(\unicode[STIX]{x1D719})$ .

Apparently, Tuck had been unaware of Tulin’s (Reference Tulin1983) work, but a chance meeting of the two occurred during a conference leading to the publication of Tuck (Reference Tuck1991b ). We are fortunate enough to possess the archived questions of the meeting, where we discover that Tulin had asked the following question:

Tuck replied that he was unsure of the generality of the reduction to problems including different geometries, but noted the connection to the Davies (Reference Davies1951) reduction (2.2):

Although Tuck’s toy reduction (1.3) should only be regarded as illustrative (the governing differential equation should rather be first order), what is apparent in his work is the desire to systematically reduce both the nonlinearity of Bernoulli’s equation, and the global effect of $\mathscr{H}$ into a single ordinary differential equation. In particular, it is Tuck’s search for a reduction of the operator, $\mathscr{H}$ , that was missing from earlier works on this topic (including Tulin Reference Tulin1983).

2.3 Other works and reductions

Certainly, Tulin and Tuck were not the only ones to seek simpler reductions or formulations of the nonlinear equations for wave–body interactions, and indeed, Tulin relates his work to the integral models proposed by Dawson (Reference Dawson1977) and Inui & Kajitani (Reference Inui and Kajitani1977). Reviews of these and other models were presented by Doctors & Dagan (Reference Doctors and Dagan1980) and Miloh & Dagan (Reference Miloh and Dagan1985), and many others. However, what distinguishes Tulin and Tuck’s work is the shifted focus towards the analytic continuation of the flow problem into the complex domain.

As we have noted in § 2.1, the low-Froude-number or low-speed limit of $\unicode[STIX]{x1D716}\rightarrow 0$ is the essential approximation in which analytical results can be derived. The subtleties of studying this singular limit can be traced back to a seminal report by Ogilvie (Reference Ogilvie1968). Ogilvie explained the importance of this regime, but then detailed certain oddities with the previously developed analytical approximations of free-surface flow past a submerged body. Chief amongst such oddities was the fact that the individual terms of a series approximation in powers of $\unicode[STIX]{x1D716}$ would fail to predict surface waves. Thus, one might decompose the surface speed into a regular series expansion and an error term, $\bar{q}$ , with

The challenge, Ogilvie realized, was that water waves are exponentially small in $\unicode[STIX]{x1D716}$ , with $\bar{q}=O(\text{e}^{-const./\unicode[STIX]{x1D716}})$ , and thus beyond all orders of any individual term of the regular series.

By linearizing about the zeroth-order solution, $q=q_{0}+\bar{q}$ for $\bar{q}\ll 1$ , and strategically preserving certain terms in Bernoulli’s equation, Ogilvie developed a general analytical approximation for the exponentially small surface waves. The approximation, however, was not asymptotically consistent, and the search for a complete numerical and analytical treatment of the low-Froude-number limit would inspire many papers in the subsequent years.

One of the key issues we will explore in this paper addresses the question of how many terms must be included in the linearization of (2.8) in order to obtain the exponential. Originally Ogilvie (Reference Ogilvie1968) had chosen $N=1$ , but later revised to $N=2$ in Ogilvie & Chen (Reference Ogilvie and Chen1982) (who quoted the work of Dagan & Tulin (Reference Dagan and Tulin1972) and in particular, the study by Keller (Reference Keller1979) in applying the WKB method to streamline ship waves).

We shall not pursue, in great detail, the history of the low-Froude-number problem that followed Ogilvie’s seminal report, but instead refer to the review papers by Tuck (Reference Tuck1991a ) and, particularly, Tulin (Reference Tulin2005) where certain aspects of the low-Froude-number difficulty are discussed. Additional historical details are presented in § 1 of Trinh, Chapman & Vanden-Broeck (Reference Trinh, Chapman and Vanden-Broeck2011), and a selection of papers on the problem is presented in table 1. The method we apply in this paper, which combines an approach of series truncation with the method of steepest descents, is unique from the previously listed works.

Note furthermore that the low-Froude-number steady-flow regime studied in this paper only forms one piece of the very large puzzle of free-surface wave–structure interactions (albeit a significant one in the context of nonlinear geometries). In other regimes, such as for the case of supercritical flows ( $\unicode[STIX]{x1D716}>1$ ), it may be possible to derive further analytical insights (cf. Craig & Sternberg Reference Craig and Sternberg1992; Tulin Reference Tulin2005; Vanden-Broeck Reference Vanden-Broeck2010). However, the connection between different regimes remains a non-trivial problem for future research; we shall return to discuss the inclusion of additional effects in § 10.

3.1 Boundary integral formulation and geometrical examples

In theory, the methods we present in this paper will apply to most general two-dimensional free-surface flows. In practice, however, in order to make analytical progress, we will be constrained by problems in which the geometry of the body is known through the specification of the angle, $\unicode[STIX]{x1D703}$ , in terms of the complex potential, $w=\unicode[STIX]{x1D719}+\text{i}\unicode[STIX]{x1D713}$ .

As specific examples, let us focus on two representative geometries: (a) flow past a localized obstruction on a channel bottom (figure 1 a) and (b) flow past a semi-infinite surface piercing ship (figure 2 a).

For flow past a varying channel bottom with dimensional upstream depth $h$ , we select the length scale $L=h/\unicode[STIX]{x03C0}$ , so that the flow region in the $w$ -plane consists of an infinite strip between $0\leqslant \unicode[STIX]{x1D713}\leqslant \unicode[STIX]{x03C0}$ . The strip is then mapped to the upper-half $\unicode[STIX]{x1D701}$ -plane, using

where under (3.4), the free surface is mapped to $\unicode[STIX]{x1D709}\geqslant 0$ , the channel bottom to $\unicode[STIX]{x1D709}\leqslant 0$ , and the flow region to $\unicode[STIX]{x1D702}\geqslant 0$ .

For flows past a semi-infinite surface piercing body, we can choose the length scale to be $L=K/U$ , where $K$ is a representative scale of the potential along the body (see (2.3) in Trinh & Chapman (Reference Trinh and Chapman2014) for further details). We assume that the free surface attaches to the body at a stagnation point, chosen without loss of generality to be $\unicode[STIX]{x1D719}=0$ . In the potential plane, the flow region consists of $\unicode[STIX]{x1D713}\leqslant 0$ . On the boundary $\unicode[STIX]{x1D713}=0$ , $\unicode[STIX]{x1D719}>0$ corresponds to the free surface and $\unicode[STIX]{x1D719}<0$ to the solid body. Since the flow is already contained within a half-plane, we do not need a further $w\mapsto \unicode[STIX]{x1D701}$ transformation, as in (3.4), but we shall set $\unicode[STIX]{x1D701}=w$ so as to use the same notation.

Bernoulli’s equation (2.1) provides an explicit relationship between the real and imaginary parts of $\unicode[STIX]{x1D6FA}$ on the free surface. However, because $\unicode[STIX]{x1D6FA}$ is analytic in the upper or lower-half $\unicode[STIX]{x1D701}$ -plane and tends to zero far upstream, there is a further Hilbert transform relationship between its real and imaginary parts. Applying Cauchy’s theorem to $\unicode[STIX]{x1D6FA}$ over a large semi-circle in the upper (channel flow) or lower (surface-piercing flow) half- $\unicode[STIX]{x1D701}$ -planes, and taking the real part gives the principal value integral,

where $j=-1$ for flow in the channel configuration and $j=1$ for flow in the ship configuration. The range of integration can be split into an evaluation over the solid boundary and over the free surface. We write

where $\mathscr{H}$ denotes the Hilbert transform operator on the semi-infinite interval,

We assume that for a given physical problem, the angle $\unicode[STIX]{x1D703}$ is known for the particular geometry, and thus the function $q_{s}$ that appears in (3.6) is known through calculating

Note that it is somewhat misleading to describe $\unicode[STIX]{x1D703}(\unicode[STIX]{x1D709})$ as ‘known’ for $\unicode[STIX]{x1D709}<0$ since in practice, we would specify $\unicode[STIX]{x1D703}$ as a function of the physical coordinates $(x,y)$ . However, specifying different forms of $\unicode[STIX]{x1D703}$ as a function of $\unicode[STIX]{x1D709}$ (or the potential, $\unicode[STIX]{x1D719}$ ) is typically sufficient to obtain the qualitatively desired geometry shape. Once the values of $q$ and $\unicode[STIX]{x1D703}$ are known on the free surface, their values everywhere within the fluid can be computed using Cauchy’s theorem. Consequently, the flow in the physical $(x,y)$ -plane can be retrieved by integrating (3.2).

For instance, we consider the step geometry,

where $0<a<b$ , which corresponds to a step of angle $\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}$ . Such topographies have been considered by King & Bloor (Reference King and Bloor1987), Chapman & Vanden-Broeck (Reference Chapman and Vanden-Broeck2006), Lustri, McCue & Binder (Reference Lustri, McCue and Binder2012), others. In this case, (3.4), (3.8) and (3.9) yields

Similarly, a semi-infinite ship with a single corner of angle $\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}$ can be specified using

Such two-dimensional hull shapes have been considered in the works of Vanden-Broeck et al. (Reference Vanden-Broeck, Schwartz and Tuck1978), Farrow & Tuck (Reference Farrow and Tuck1995), Trinh et al. (Reference Trinh, Chapman and Vanden-Broeck2011), and others. Choosing the dimensional length, $L=K/U$ where $K$ is the value of the potential at the corner sets its non-dimensional position to $\unicode[STIX]{x1D719}=-1$ . Then using (3.8) and (3.11), the $q_{s}$ function is given by

Recall the ship problem does not require an additional mapping to the half plane, so we can set $\unicode[STIX]{x1D701}=w$ to preserve the notation.

In this paper, we will only consider geometries which contain strong singularities in the complex plane; these correspond to poles or branch points in $q_{s}$ . For instance, the step (3.10) and ship (3.12) contain singularities at the corners and stagnation points of the body. We will see in § 8 that these singularities are often responsible for the creation of the surface waves. Weaker singularities, such as the case of flow past a smoothed hull with a jump in the curvature, will be the subject of a forthcoming paper (cf. also § 7 of Trinh & Chapman Reference Trinh and Chapman2014). Flows past closed objects also present challenging cases for study because typically the geometry is not known as a function of the potential variables, and there are further difficulties with the prediction of wake separation points (figure 1 b).

The situation of steady-state gravity waves requires specification of additional radiation or boundary conditions. In the case of step flow in a channel, we impose a uniform flow upstream,

For the case of stern flow past a ship, we require that the flow separates at a stagnation point, and the necessary conditions on $q$ and $\unicode[STIX]{x1D703}$ were derived by Dagan & Tulin (Reference Dagan and Tulin1972). If it is the case that $\unicode[STIX]{x1D70E}\geqslant 1/3$ , then the flow separates horizontally from the body and this requires that

If $\unicode[STIX]{x1D70E}<1/3$ , then the flow separates so as to make an interior angle of $2\unicode[STIX]{x03C0}/3$ with the body. Steady-state bow flows are typically not possible unless the free surface is allowed to form a splash. Further extended discussions of the requisite boundary conditions are found in Vanden-Broeck (Reference Vanden-Broeck2010), Trinh et al. (Reference Trinh, Chapman and Vanden-Broeck2011), and the references therein.

Thus the governing equations consist of Bernoulli’s equation (2.1), the boundary integral equation (3.6), a specified geometry, and the radiation or boundary conditions (3.13).

4.1 Divergence of the asymptotic approximation

The leading-order approximation, $q_{0}$ , contains singularities in its analytic continuation into the complex plane. For instance, in the case of the step geometry (3.10), we see that $q_{0}$ contains branch points at $\unicode[STIX]{x1D709}=-a$ and $\unicode[STIX]{x1D709}=-b$ , which correspond to the corner and stagnation points.

Because of the singularly perturbed nature of Bernoulli’s equation in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ , we can see that at each order in the asymptotic procedure, the calculation of $q_{n}$ and $\unicode[STIX]{x1D703}_{n}$ will depend on differentiation of the previous order, $q_{n-1}$ and $\unicode[STIX]{x1D703}_{n-1}$ . Since the leading-order approximation is singular, this has the effect of increasing the power of the singularities with each subsequent term. Thus in the limit that $n\rightarrow \infty$ , $q_{n}$ and $\unicode[STIX]{x1D703}_{n}$ will diverge. As argued in the works of e.g. Chapman, King & Adams (Reference Chapman, King and Adams1998), Chapman & Vanden-Broeck (Reference Chapman and Vanden-Broeck2006), the divergence at $O(\unicode[STIX]{x1D716}^{n})$ is captured through a factorial over power ansatz,

where $\unicode[STIX]{x1D712}=0$ at the particular singularity. The functions $Q$ , $\unicode[STIX]{x1D6E9}$ , and $\unicode[STIX]{x1D712}$ may be regarded as either functions of $w$ or of $\unicode[STIX]{x1D701}$ . If there are multiple singularities contributing to the divergence (see Trinh & Chapman (Reference Trinh and Chapman2014) for an example), then we must include a summation over similar ansatzes of the form (4.2) where $\unicode[STIX]{x1D712}=0$ at each individual singularity. There exist other cases where a more general form of the divergence is required, and this is documented in Trinh & Chapman (Reference Trinh and Chapman2015).

In order to derive the components, $Q$ , $\unicode[STIX]{x1D712}$ , and $\unicode[STIX]{x1D6FE}$ , we can examine the $O(\unicode[STIX]{x1D716}^{n})$ form of the system, and take the limit $n\rightarrow \infty$ . In this paper, the late-orders behaviour is not a crucial part of the analysis, but we collect the functional forms of the components in appendix A.

5.1 Further notes on the analytic continuation process

There are subtle aspects of analytically continuing $q$ and $\unicode[STIX]{x1D703}$ off the free surface, and the relationship of these quantities with the usual physical complex velocity given by $\operatorname{d}w/\operatorname{d}z=u-\text{i}v=q\text{e}^{-\text{i}\unicode[STIX]{x1D703}}$ . The $q_{c}$ and $\unicode[STIX]{x1D703}_{c}$ that appear on the left-hand side of (5.5) are complex valued, but reduce to the physical speed and velocities on the free surface, $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}>0$ . Although their combined values, $q_{c}\text{e}^{-\text{i}\unicode[STIX]{x1D703}_{c}}$ , are related to the complex velocity, $u-\text{i}v$ , their individual values are not known without further work.

As an example of this particular subtlety, we may consider the analytic continuation of $\unicode[STIX]{x1D703}$ , evaluated along the physical boundary, $\unicode[STIX]{x1D701}\leqslant 0$ . In general, this value, $\unicode[STIX]{x1D703}_{c}$ , will not be the physical angle as defined in (3.9) or (3.11). Instead, the analytic continuations of the individual components, $q_{c}$ and $\unicode[STIX]{x1D703}_{c}$ , are related to the physical angle through $\operatorname{Im}[\log (q_{c}\text{e}^{-\text{i}\unicode[STIX]{x1D703}_{c}})]$ .

In this work, we will depend on asymptotic derivations of the analytic continuation variables, $q_{c}$ and $\unicode[STIX]{x1D703}_{c}$ . However, it is also possible to numerically solve the nonlinear equations (5.6) for such quantities; this is done for the study of travelling waves of permanent form in Crew & Trinh (Reference Crew and Trinh2016).

7.1 Truncation and optimal truncation

The crucial decision in studying the system (7.2) is to choose how many terms to include in the truncated series, $q_{r}$ and $\unicode[STIX]{x1D703}_{r}$ , which determines the forcing function $R$ . That is, we must choose the truncation value of $N$ .

In early studies of the low-Froude-number problem, researchers had used $N=1$ , and then later corrected to $N=2$ . However, at a fixed value of $\unicode[STIX]{x1D716}$ , the truncation error in the divergent asymptotic expansion will continue to decrease for increasing values of $N$ . It reaches a minimum at an optimal truncation point, $N={\mathcal{N}}$ , and diverges afterwards. For small $\unicode[STIX]{x1D716}$ , the optimal truncation point is found where adjacent terms in the series are approximately equal in size, or $|\unicode[STIX]{x1D716}^{{\mathcal{N}}}q_{{\mathcal{N}}}|\sim |\unicode[STIX]{x1D716}^{{\mathcal{N}}-1}q_{{\mathcal{N}}-1}|$ . Using the divergent form (4.2), we find that the optimal truncation point is given by

Thus, for a fixed point in the domain, $w$ , and in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ , the optimal truncation point, ${\mathcal{N}}$ , tends to infinity, and we must include infinitely many terms of the series. In § 9, we will examine the effect of different truncations on the comparison between numerical solutions and asymptotic approximations.

Most emphatically, we remark that apart from the hidden $O(\unicode[STIX]{x1D716}^{2})$ terms multiplying $\bar{q}$ in the brackets of (7.2a ), this reduced equation is an exact result to all orders in $\unicode[STIX]{x1D716}$ . We have only employed the fact $\bar{q},\bar{\unicode[STIX]{x1D703}}\ll 1$ , but the inhomogeneous term, $R(w)$ is exact. When the regular series expansions are optimally truncated, $\bar{q}$ and $\bar{\unicode[STIX]{x1D703}}$ are exponentially small, and consequently, the $O(\bar{q}^{2},\bar{\unicode[STIX]{x1D703}}^{2})$ terms are not needed to derive the leading-order behaviour of the exponentials.

Though it is now written as a first-order differential equation, the system (7.2) is difficult to solve, since it involves real and complex components, in addition to the integral transforms embedded in the inhomogeneous term. The key in the following sections will be to argue that the Hilbert transform can be neglected.

8.1 Integral form of the solution

To begin, we integrate (7.2b ) to obtain

In the above integral, the initial point of integration, $w_{0}$ , is chosen to be the particular point that causes the divergence of the late-order terms in (4.2). Thus, $\unicode[STIX]{x1D712}$ defined in (7.2) matches the definition of the singulant function in the late-orders ansatz.

Next integrating (7.2c ),

where $-\log \unicode[STIX]{x1D6EC}$ is the constant of integration, and $w^{\ast }$ can be chosen wherever the integral is defined. We also note that $P_{1}$ given by (8.2) is related to $Q$ given in (A 3). We can then write

Solving (7.2a ) now yields a solution for the integro-differential equation (though note that it is not quite explicit due to the reliance on the Hilbert transform),

where we have introduced the integral

which will be used to extract the prefactor of the exponential. The start point, $w_{s}$ , can be chosen based on additional boundary or radiation conditions. For channel flow, the water wave problem will naturally impose a radiation condition which requires the free surface to be flat at $\unicode[STIX]{x1D719}=-\infty$ (upstream), so taking $w_{s}\rightarrow -\infty$ , the constant of integration in (8.4a ) is zero. Similarly, for the surface-piercing problem, $w_{s}$ can be taken to be the stagnation point, $w_{s}=0$ , where $q=0$ .

8.2 The steepest descent paths

For the integral (8.4b ), the paths of steepest descent are given by constant contours of $\operatorname{Im}\unicode[STIX]{x1D712}$ . Let us envision a generic situation in which the steepest descent topology is as shown in the figure 4(a). A critical point (branch point, singularity, or saddle point) of the integrand is assumed to lie at $w=w_{0}$ . The valleys, where $\text{Re}[\unicode[STIX]{x1D712}(w)]\leqslant \text{Re}[\unicode[STIX]{x1D712}(w_{0})]$ , are shown shaded. In the illustrated case, the path of steepest descent from the initial, $t=w_{s}$ , and final point, $t=w$ , is joined by a single valley of the integrand. Thus, we write

and the relevant contributions to the integral solely depend upon an expansion about the endpoints. When we perform this analytically in § 8.3, we find that the $I_{endpoint}$ contribution yields further algebraic orders of the base asymptotic series, $q_{r}$ and $\unicode[STIX]{x1D703}_{r}$ . In other words, (8.5) leads to further deflections of the free surface due to the moving body, but does not produce an oscillatory component.

Figure 4.

Illustration of steepest descent contours for a problem with a single wave-generating singularity, $w=w_{0}$ . (a) When the endpoint of integration lies to the left of the Stokes line, $I\sim I_{endpoints}$ ; (b) when the end of integration crosses the Stokes line, the steepest descent path must include the contribution from the saddle.

The oscillatory contributions arise from the deformation of the path of integration through the critical points of the integrand. A typical situation is illustrated in figure 4(b). There, we see that if the endpoint, $t=w$ , crosses the Stokes line, $\operatorname{Im}[\unicode[STIX]{x1D712}]=\operatorname{Im}[\unicode[STIX]{x1D712}(w_{0})]$ from the point, $t=w_{0}$ , then the steepest descent paths from each individual endpoint tend to distinct valleys. These valleys must be joined by a contour through the saddle or around the critical point. In such cases, we shall write

In this work, we do not discuss the global properties of the steepest descent curves and the associated Stokes lines. These results require careful consideration of the Riemann-sheet structure of the integrand, and can be found in the accompanying work (Trinh Reference Trinh2016). Instead, we assume that the integral is approximated by (8.6), and numerically verify the approximation in § 9.

8.3 Endpoint contributions

It can be verified that the dominant contributions from the endpoints will re-expand the regular perturbation series to higher orders. In particular, we see from inspection that locally near $t=w$ , the integrand is exponentially large and of order $\text{e}^{\unicode[STIX]{x1D712}/\unicode[STIX]{x1D716}}$ , and this serves to cancel the exponentially small factor in (8.4a ).

The form of the integral in (8.4) is unwieldy, so we shall demonstrate this re-expansion for the simplest case where the regular series, $q_{r}$ and $\unicode[STIX]{x1D703}_{r}$ , contains only a single term. Setting $q_{r}=q_{0}=q_{s}$ and $\unicode[STIX]{x1D703}_{r}=\unicode[STIX]{x1D703}_{0}=0$ from (4.1a,b ) and using (7.2d–f ) the integral is reduced to

Integrating by parts, we have

where we have assumed that the boundary term contributions due to $\unicode[STIX]{x1D711}=w_{s}$ are zero from the boundary or radiation conditions (see the discussion following (8.4b )). Repeated integration by parts to the integral term of (8.7) will extract further contributions. We keep only the first term above, and recombine with $\bar{q}$ in (8.4a ) to obtain

where we have used (5.3) and (7.2b ). Matching real and complex parts then yields $\bar{q}=jq_{0}\mathscr{H}[\bar{\unicode[STIX]{x1D703}}]$ and $\bar{\unicode[STIX]{x1D703}}=-\unicode[STIX]{x1D716}q_{0}^{2}q_{0}^{\prime }$ , which are the first correction terms derived in (4.1c,d ). The procedure to obtain further corrections proceeds in a similar way using repeated integration by parts, and thus producing further powers of $\unicode[STIX]{x1D716}$ .

Including more terms from $q_{r}$ and $\unicode[STIX]{x1D703}_{r}$ and performing the re-expansion procedure about the endpoints will follow the same idea; algebraically, however, it is much easier to work with the integro-differential equation (7.2a ), or simply with the original system of two coupled equations, as was done in § 4.

8.4 The Hankel contour

In the limit $\unicode[STIX]{x1D716}\rightarrow 0$ , the contributions to the contour integral associated with wave motion arises from cases where the contour is deformed past critical points in the integrand, or near saddle points, where $\unicode[STIX]{x1D712}^{\prime }=0$ . Let us assume that the contribution occurs at a point, $w=w_{0}$ , where $\unicode[STIX]{x1D712}^{\prime }(w_{0})=0$ and that $w_{0}$ lies off the free surface (see § 10.2 for comments on other cases). Locally near this point, we assume that the leading-order speed behaves as

for some constants $c$ and $\unicode[STIX]{x1D6FC}$ . Consequently, it follows from (7.2b ) that

One of the key issues that underlies the Tulin reduction of the water-wave equations concerns the role of the non-local Hilbert transform, $\hat{\mathscr{H}}$ , which appears in the forcing function $R(w;\hat{\mathscr{H}}[\bar{\unicode[STIX]{x1D703}}])$ in (7.2d ). In the previous section, we discovered that the Hilbert transform plays an important role in the further development of the endpoint contributions. Indeed, it appears as a contributing term in each order of the asymptotic process, beginning from $q_{1}$ in (4.1c,d ).

Upon applying the method of steepest descents to the integral $I$ , we had separated the contributions into those due to endpoints and those due to saddle points (assumed to lie away from the free surface). Let us assume that the contribution due to the saddle point is exponentially small along the free surface, with $\bar{q}_{exp}=O(\text{e}^{-\unicode[STIX]{x1D712}/\unicode[STIX]{x1D716}})$ in (8.4a ), with $\text{Re}(\unicode[STIX]{x1D712})>0$ on the free surface, and similar relations for $\bar{\unicode[STIX]{x1D703}}$ . As the variable of integration approaches the saddle point, $\unicode[STIX]{x1D711}\rightarrow w_{0}$ , the term

that appears in $R$ of (7.2d ) involves the integration of an exponentially small term along the free surface, where it remains exponentially small. Thus, this term is subdominant to the square-bracketed contributions in (7.2d ). This argument does not work for the endpoint contributions; the Hilbert transform is not negligible, as it involves the integration of the regular perturbative series via $\hat{\mathscr{H}}[\unicode[STIX]{x1D703}_{r}]$ .

We now approximate the contribution near the saddle points. As discussed in § 7.1, as $\unicode[STIX]{x1D716}\rightarrow 0$ , the optimal truncation point, $N\rightarrow \infty$ . Thus using the divergent form (4.2), we have that

away from the real axis. Recall that the components, $Q$ , $\unicode[STIX]{x1D6FE}$ , and $\unicode[STIX]{x1D712}$ are given in appendix A. Thus substituting (8.13) into (8.4b ), we have

In (8.14), the contour $C_{0}$ corresponds to the steepest descents contour in the vicinity of the critical point, $\unicode[STIX]{x1D711}=w_{0}$ . Near this point, we simplify the exponential by making a coordinate transformation with

from (8.11). It thus follows from (8.14) that

Based on the steepest descent topologies corresponding to the ship and step problems, the relevant local contribution near $\unicode[STIX]{x1D711}=w_{0}$ occurs in the form of a Hankel contour – that is, an integral about the branch cut $u\in \mathbb{R}^{-}$ , beginning from $u=-\infty -0\text{i}$ , looping around the origin, and tending to $u=-\infty +0\text{i}$ . These steepest descent topologies are the subject of the companion paper (Trinh Reference Trinh2016). By the integral definition of the gamma function,

Therefore, combining with (8.16), we have $I_{exp}\sim 2\unicode[STIX]{x03C0}\text{i}\unicode[STIX]{x1D716}^{1-\unicode[STIX]{x1D6FE}}$ . We may now substitute this approximation for $I_{exp}$ into (8.4a ) to obtain the part of $\bar{q}$ that is switched-on by the saddle point contribution. This yields

where $Q(w)$ is given in (8.3).

Recall that the functional form of integrand in (8.4b ) changes in the limit the real axis is approached due to the presence of the complex Hilbert transform in $R$ . Moreover, the derivation of (8.18) only takes in account the steepesst descent paths for analytic continuation into the upper half- $\unicode[STIX]{x1D701}$ -plane ( $k=1$ ). Repeating the process for the lower half reflects the steepest descent paths about the real axis.

Rather than repeating the procedure and observing the change in signs of $k$ , we can note that wave component of $q$ or $\unicode[STIX]{x1D703}$ must be real along the physical free surface. Consequently, by the Schwarz reflection principle, the leading-order wave solution must be the sum of (8.18) and its complex conjugate. Thus on the free surface, where $w=\unicode[STIX]{x1D719}$ ,

In numerical computations, it is typically easiest to compare the wave amplitudes in the far field, where $w=\unicode[STIX]{x1D719}\rightarrow \infty$ , and the mean speed, $q_{r}\rightarrow 1$ . For simplicity we choose the integral limit $w^{\ast }$ in $Q$ to be on the free surface. Substituting (8.3) into (8.19) and simplifying gives the amplitude of the oscillations,

where the constant $\unicode[STIX]{x1D6EC}$ appears in the expression for $Q$ and its value must be computed numerically for general nonlinear problems (see appendix A). In the end, (8.20) provides a closed-form formula for the amplitude of the waves, which only depends on specification of the two solutions, $q_{0}$ and $q_{1}$ , calculated from (4.1), and late-order components, $\unicode[STIX]{x1D6EC}$ , $\unicode[STIX]{x1D6FE}$ , and $\unicode[STIX]{x1D712}$ , given in appendix A. Note that $q_{1}$ is needed to calculate the value of $\unicode[STIX]{x1D6EC}$ .

Finally, we can establish an equation between the exponentially small waves in $q$ and $\unicode[STIX]{x1D703}$ . Upon substituting (7.1) into (5.6b ), we have that

9.1 Truncation at other values of $N$

The result of the Hankel integration, in (8.19), which depends on the late-terms behaviour of the ansatz (4.2), was derived on the assumption that truncation of the regular series expansion occurs at optimal truncation with $N\rightarrow \infty$ as $\unicode[STIX]{x1D716}\rightarrow 0$ . In practice, however, it may be useful to develop reduced models that truncate at pragmatic values of $N$ (i.e. one or two); let us discuss what loss in accuracy results from this.

As we have shown, the surface speed is expressed as

where we have written the exponential (8.18) in a form so as to separate a numerical prefactor, ${\mathcal{A}}$ , from the functional dependence, $F(\unicode[STIX]{x1D719})$ .

The case of truncation at $N=1$ is asymptotically inconsistent since, if $q=q_{0}+\bar{q}$ , then $\bar{q}$ contains both an $O(\unicode[STIX]{x1D716})$ error and also the exponentially small waves we desire. But from (7.2a ) and (7.2c ), we see that it would then be incorrect to ignore $O(\bar{q}^{2})$ errors in deriving the linearized equation. Despite these issues it is still instructive to examine the truncation at $N=1$ . This yields

and the solution

for an initial point of integration, $t=w_{s}$ , where $\bar{q}(w_{s})\rightarrow 0$ .

Thus the exponential argument, $\text{e}^{-\unicode[STIX]{x1D712}/\unicode[STIX]{x1D716}}$ , of the solution can be derived from (9.2), and indeed the steepest descent paths of § 8 are still applicable. However we have inaccurately predicted the functional dependence, $F(\unicode[STIX]{x1D719})$ .

The case of truncation at $N=2$ is much more interesting. It yields

Thus, comparing (9.4) to (7.2), we have obtained the precise left-hand side required. In developing the leading-order exponential, the only error in (9.4) is due to the replacement of the exact right-hand side by its $N=2$ truncation. This only affects the integrand of (8.4), whose role was to determine the prefactor ${\mathcal{A}}$ . In other words, the linear equation (9.4) allows us to develop nearly all of the leading-order exponential. We will see in the numerical computations that this equation provides a very close fit to results from the full nonlinear model. The two-term truncation is the closest analogy to the Tulin (Reference Tulin1983) model (see § 10 for an extended discussion), as well as various other models (cf. the review by Doctors & Dagan Reference Doctors and Dagan1980). However the fact that such two-term approximations do not accurately predict the prefactor ${\mathcal{A}}$ is not well established in the literature.

9.2 A simplified nonlinear model for $N=2$

An even simpler formulation to (9.4) can be developed at the risk of slightly increased inaccuracy. This form serves as an extremely useful toy model due to the complete removal of the Hilbert transform, and thus in this sense, it provides the strongest analogy to the Tuck reduction of (2.6). Let us return to combined integro-differential equation in (5.8). Approximating $\text{e}^{\;j\hat{\mathscr{H}}[\unicode[STIX]{x1D703}]}\sim 1$ , we obtain

In order to avoid confusion, we shall write the solution of (9.5) as $\tilde{q}$ . Expressing $\tilde{q}$ as a regular series expansion $\tilde{q}=\tilde{q}_{0}+\unicode[STIX]{x1D716}\tilde{q}_{1}+\cdots \,$ , we find

The leading-order solution is as expected, but ignoring the $\hat{\mathscr{H}}[\unicode[STIX]{x1D703}]$ complex Hilbert transform has the effect of changing the correction term, $\tilde{q}_{1}$ . If we return to the correct $q_{1}$ in (4.1c ), we find that

upon using (5.3) and the solution for $\unicode[STIX]{x1D703}_{1}$ in (4.1d ). Recall that near the wave-generating singularity, $\hat{\mathscr{H}}[\unicode[STIX]{x1D703}_{1}]$ remains bounded, but $q_{0}$ and its derivatives are singular according to (8.10). In other words, the simplified formulation of (9.5) has replaced $q_{1}$ in the full model with its local behaviour near the singularity.

Substituting $q=\tilde{q}_{0}+\unicode[STIX]{x1D716}\tilde{q_{1}}+\cdots +\bar{q}$ into the simplified nonlinear model (9.5) and linearizing in $\bar{q}$ yields

where we have withheld the right-hand side, $\tilde{R}$ , for clarity. From (8.2) and (9.7), we can verify that as $w\rightarrow w_{0}$ , $P_{1}^{\prime }\sim 5q_{0}^{\prime }/q_{0}$ . Comparing the bracketed terms in (9.4) and (9.8), we thus see that in the simplified model, the $\hat{\mathscr{H}}[\unicode[STIX]{x1D703}]$ terms are neglected but we have nevertheless preserved the inner limit of (9.4). Since the Hankel contour analysis of § 8.4 depends only on local properties near the relevant singularities, we conclude that the salient details of the water wave analysis are still largely preserved for the model in (9.5). Consequently, as it concerns the resultant waves, the model is not able to duplicate the exact functional form of the prefactor, but does produce a prefactor that contains the correct inner limit.

In the study of Trinh & Chapman (Reference Trinh and Chapman2015), this simplified nonlinear problem was used as a toy model for the study of wave–structure interactions with coalescing singularities. By duplicating the derivation of § 8, the analogous formula to (8.20) can be developed for (9.5).

A detailed summary of the full and simplified models we have presented thus far is shown in table 2.

Table 2.

A comparison of different truncated models for the study of gravity waves past moving bodies, and their ability to correctly predict (yes/no) the leading-order exponential, ${\mathcal{A}}F(\unicode[STIX]{x1D719})\text{e}^{-\unicode[STIX]{x1D712}/\unicode[STIX]{x1D716}}$ (to be added to its complex conjugate on the free surface). The case of $N=2^{\ast }$ is a two-term truncation that extracts the singular behaviour of the Hilbert transform, and thus obtains the correct functional form of $F(\unicode[STIX]{x1D719})$ near the relevant singularity, $w=w_{0}$ .

9.3 Numerical results of the stern-flow problem

We now compare numerical solutions of three reduced models with the full nonlinear stern-wave problem in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ . The models include (i) the full nonlinear problem (5.8), or more conventionally, the solution of Bernoulli’s equation (2.1) and the boundary integral equation (3.6); (ii) the $N=2$ truncated linear model in (9.4); and (iii) the $N=2^{\ast }$ nonlinear model in (9.5).

We study the case of a semi-infinite rectangular ship given by (3.11) and (3.12) with $\unicode[STIX]{x1D70E}=1/2$ , and thus

which is the leading-order speed associated with a hull with a right-angled corner ( $3\unicode[STIX]{x03C0}/2$ in the interior of the fluid) at $w=-1$ , and a stagnation point at $w=0$ .

The solution is computed in each of the three cases, and the amplitude of the water waves far downstream, with $w=\unicode[STIX]{x1D719}\gg 1$ , is extracted. For the two models with truncations at $N=2$ and $N=2^{\ast }$ , recall that the real-valued solution on the axis is formed by adding the complex-valued solution to its complex conjugate [see (8.19) and the surrounding discussion]. Thus, for these two cases, the amplitude of $q$ is taken to be twice the amplitude of $\text{Re}(q_{exp})$ .

In order to solve the full water-wave equations (5.8), we use the numerical algorithm described in Trinh et al. (Reference Trinh, Chapman and Vanden-Broeck2011). In brief, a stretched grid is applied near the stagnation point, and a finite-difference approximation of the boundary integral is calculated using the trapezoid rule. At a singular point of the integral, a quadratic interpolant is applied between the point and its two neighbours, and the resultant quadratic is calculated exactly. For more details of the numerical scheme see, e.g. Chapter 7 of Vanden-Broeck (Reference Vanden-Broeck2010) and the references therein. From the predicted wave amplitude (8.20), we have

where $\tilde{C}\approx 2.215$ , $\unicode[STIX]{x1D6FE}=6\unicode[STIX]{x1D70E}/(1+3\unicode[STIX]{x1D70E})$ , and here, $\unicode[STIX]{x1D70E}=1/2$ . The numerical prefactor $\tilde{C}$ requires a generic calculation (cf. figure 10 of Trinh et al. Reference Trinh, Chapman and Vanden-Broeck2011). Both numerical amplitude measurements (stars), as well as the asymptotic prediction (upper dashed line), are shown in figure 5.

The $N=2$ and $N=2^{\ast }$ truncated models can be solved as initial value problems. Due to the singular nature of the stagnation point, where $q=O(w^{1/2})$ , we use the coordinate transformation $s(w)=w^{1/2}$ , and solve the associated differential equations in $s$ . The asymptotic behaviour (9.9) is used to provide the initial value for $q$ at a point near $s=0$ . In the simulations, we typically used $s=10^{-5}$ for this point, and the resultant amplitudes are verified to be independent of the initial condition.

The combined results are shown in figure 5. The leading-order asymptotic approximations fit the data closely, and both the simplified nonlinear ( $N=2^{\ast }$ , shown as circle) and linear ( $N=2$ , shown as stars) formulations duplicate the requisite behaviours reviewed in table 2. Note that although the emphasis of this work is largely on the development of the mathematical framework, the numerical results indicate that the low-speed theory is able to predict both the exponential dependence to some of the typical values of the Froude number in ship hydrodynamical applications (see Wehausen (Reference Wehausen1973) for more details).

Figure 5.

Comparison of wave amplitudes for the full nonlinear model, the $N=2$ truncated linear model, and the $N=2^{\ast }$ truncated nonlinear model. For the two truncated models, the amplitude is multiplied by two to account for the analytic continuation. The dashed lines correspond to leading-order asymptotic approximations of the full nonlinear model (top) and the simplified nonlinear model (bottom).

9.4 Numerical results of the step-wave problem

Numerical solutions for flow past a step of angle $\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}$ are computed by solving the coupled system of (2.1) and (3.6), but now with the adjusted map (3.4) and leading-order speed (3.10). Description of the numerical scheme and asymptotic comparisons with the optimally truncated model have been presented in Chapman & Vanden-Broeck (Reference Chapman and Vanden-Broeck2006, figures 5–8) and Lustri et al. (Reference Lustri, McCue and Binder2012, figures 4, 9). In this section, we provide some additional discussion of the Stokes line structure, and how it relates to the formation of waves.

A typical solution is shown in figure 6(a) for the case of a rectangular step with $\unicode[STIX]{x1D70E}=1/2$ at $\unicode[STIX]{x1D716}=0.35$ . Note the sharpened crests and wide troughs of the surface waves; we have chosen the Froude number to be particularly large in order to better visualize the free surface. In order to plot streamlines in the physical plane, we calculate values of the complex velocity, $q\text{e}^{-\text{i}\unicode[STIX]{x1D703}}$ , in the lower half- $\unicode[STIX]{x1D701}$ -plane using the analytic continuation of the boundary integral (5.5) with $j=-1$ and $k=-1$ . Once $q$ and $\unicode[STIX]{x1D703}$ are known, the physical $x$ and $y$ -coordinates along streamlines follow from integrating (3.2).

Figure 6.

(a) Flow past the rectangular step defined by (3.10) with $\unicode[STIX]{x1D70E}=1/2$ , $a=1$ and $b=2.2$ . The solution is computed using the full nonlinear equations (2.1) and (3.6) at $\unicode[STIX]{x1D716}=0.35$ . Streamlines are calculated by analytic continuation of the boundary integral into the fluid region. Also shown is the projection of the Stokes line, where $\operatorname{Im}\unicode[STIX]{x1D712}=0$ and $\text{Re}\unicode[STIX]{x1D712}\geqslant 0$ , which emerges from the corner, $\unicode[STIX]{x1D701}=-a$ or $(x,y)\approx (0,0.86)$ . An enlargement of the region near the Stokes line is shown in (b). The drag force on the step is shown in (c) and computed from (9.11) for different $\unicode[STIX]{x1D716}$ . As $\unicode[STIX]{x1D716}\rightarrow 0$ , the drag behaves as $D=-\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}(\sqrt{a}-\sqrt{b})^{2}/\sqrt{4ab}+O(\unicode[STIX]{x1D716}^{2})$ .

In § 8.2, we showed that for a given flow problem, Stokes lines within the integration plane of (8.4b ) are associated with sudden changes in the steepest descent paths. Across such lines, waves are switched-on via a Hankel contour around a particular singularity (cf. § 8.4). We also noted in § 5.1 that, although Stokes lines lie in a non-physical region of space (the analytic continuation of the free surface), they can often be projected into physical space. In the case of the rectangular step of figure 6(a), we do this computation as follows.

First, a separate analysis indicates that the stagnation point of the step does not produce a Stokes line (Chapman & Vanden-Broeck Reference Chapman and Vanden-Broeck2006). Thus, the initial point of integration of $\unicode[STIX]{x1D712}$ in (8.1) is chosen at the corner of the step, $\unicode[STIX]{x1D701}=-a$ or $w_{0}=-\text{log}|a|-\unicode[STIX]{x03C0}\text{i}$ . The Stokes line then consists of points where $\{\operatorname{Im}\unicode[STIX]{x1D712}=0,\text{Re}\unicode[STIX]{x1D712}\geqslant 0\}$ , and this is computed by numerical quadrature in the $\unicode[STIX]{x1D701}$ -plane. Once the Stokes line is known in the $\unicode[STIX]{x1D701}$ -plane or equivalently the $w$ -plane, it can be projected into the $(x,y)$ -plane by interpolation of the streamlines of the flow. The result is shown in figure 6(a) and a close-up in figure 6(b). The Stokes line marks the location in the flow where the waves switch on following the rise in the step. At finite values of $\unicode[STIX]{x1D716}$ , the transition occurs in a boundary layer of appreciable thickness (Berry Reference Berry1989).

There are several interesting consequences of considering the low-speed limit of $\unicode[STIX]{x1D716}\rightarrow 0$ and viewing the Stokes line projected into the physical plane. For example, King & Bloor (Reference King and Bloor1987) had considered the non-dimensional drag force on the face of the step given by,

where $h$ is the step height in the physical plane. In the limit $\unicode[STIX]{x1D716}\rightarrow 0$ , the drag can be approximated by integrating (9.11) with the speed given by the leading-order speed, $q_{0}$ , and this is verified in figure 6(c). The higher-order terms, $q_{1}$ , $q_{2}$ , and so forth produce further algebraic corrections of the drag force. In the case of general angled steps, waves may be switched on in the region between $\unicode[STIX]{x1D701}=-b$ and $\unicode[STIX]{x1D701}=-a$ , and therefore, the wave field may affect the drag force to an exponentially small degree (cf. Stokes line diagrams in Lustri et al. (Reference Lustri, McCue and Binder2012) and Trinh (Reference Trinh2016)).