Gravity-capillary waves in reduced models for wave-structure interactions

TL;DR

This paper introduces a reduced model with a Sommerfeld-type boundary condition to accurately compute steady gravity-capillary waves in low-speed flows, addressing longstanding numerical challenges in wave-structure interaction simulations.

Contribution

It presents a novel asymptotic reduction and boundary condition approach for steady gravity-capillary waves, enabling more accurate numerical solutions in low-speed regimes.

Findings

Validated low-speed wave predictions using exponential asymptotics

Demonstrated sensitivity of full equations to boundary condition errors

Provided a practical method for resolving the radiation problem

Abstract

In order to determine the steady-state subcritical gravity-capillary waves that are produced by potential flow past a wave-making body, it is typically necessary to impose a radiation condition that allows for capillary waves upstream, but disallows those corresponding to gravity. However, this radiation condition is not known a priori and consequently, the computation of accurate numerical solutions to the steady-state problem remains problematic. Although the physical model might be modified (e.g. with viscosity), recovery of the original problem is not always possible. In this work, we discuss the above radiation problem, and show how, in the low-speed limit, the steady gravity-capillary waves can be resolved using a Sommerfeld-type boundary condition applied to an asymptotically reduced set of water-wave equations. These results allow us to validate the specialized classes of low-speed waves theoretically predicted by Trinh & Chapman (2013) using methods in exponential asymptotics [J. Fluid Mech. 724, pp. 392--424]. The issues of numerically solving the full set of gravity-capillary equations for a potential flow are discussed, and the sensitivity to errors in the boundary conditions is clearly demonstrated.

Equations

$$\mathrm{F}=\frac{U}{\sqrt{gL}}\text{and}\mathrm{T}=\frac{\sigma }{\rho g{L}^2}.$$

$$U^2=\left[\frac{g}{\tilde{k}}+\frac{\sigma }{\rho }\tilde{k}\right]\tanh (\tilde{k}h).$$

$${\mathrm{F}}^2=\left[\frac{1}{k}+\mathrm{T}k\right]\tanh (k)\equiv G(k;\mathrm{T}).$$

$$k=\tilde{k}{h}_{\text{up}}.$$

$$h_{\text{d/u}}=\frac{h_{\text{down}}}{h_{\text{up}}}\text{and}{U}_{\text{d/u}}=\frac{U_{\text{down}}}{U_{\text{up}}},$$

$$\mathrm{F}={\mathrm{F}}_{\text{down}}\frac{\sqrt{h_{\text{d/u}}}}{U_{\text{d/u}}},\mathrm{T}={\mathrm{T}}_{\text{down}}{h}_{\text{d/u}}^2,k=\frac{k_{\text{down}}}{h_{\text{d/u}}},$$

$${\mathrm{F}}^2=h_{\text{d/u}}^2\left[\frac{1}{k}+\mathrm{T}k\right]\tanh (h_{\text{d/u}}k)\equiv H(k;\mathrm{T},h_{\text{d/u}}).$$

$${\mathrm{F}}^2=\beta ϵ\text{and}\mathrm{T}=\beta \tau {ϵ}^2,$$

$$U_{\text{d/u}}=\sqrt{b}+\mathcal{O}(ϵ)\sim \frac{1}{h_{\text{d/u}}}.$$

$$\beta \tau {ϵ}^2{k}^2-\beta ϵk+1=0$$

$$k_{\text{up}}=\frac{1\pm \sqrt{1-A}}{2\tau ϵ},$$

$$k_{\text{down}}=k_{\text{up}}{h}_{\text{d/u}}=\frac{b\pm \sqrt{b^2-A}}{2\tau \sqrt{b}ϵ}.$$

$$:$$

$$L=\frac{h_{\text{up}}}{\pi },$$

$$\zeta =\xi +\mathrm{i}\eta ={\mathrm{e}}^{-w}={\mathrm{e}}^{-(\varphi +\mathrm{i}\psi )}.$$

$$\frac{dw}{dz}=q{\mathrm{e}}^{-\mathrm{i}\theta },$$

$$\log q(\xi )=-\frac{1}{\pi }\left[{\int }_{-\infty }^0\frac{\theta ({\xi }^{\prime })}{{\xi }^{\prime }-\xi }d{\xi }^{\prime }+\text{vbox}-{\int }_0^{\infty }\frac{\theta ({\xi }^{\prime })}{{\xi }^{\prime }-\xi }d{\xi }^{\prime }\right].$$

$$\theta (\xi )=\{\begin{array}{ll}0& -\infty <\xi <-b,\\ \pi /2& -b<\xi <-1,\\ 0& -1<\xi <0,\end{array}$$

$$-\frac{1}{\pi }{\int }_{-\infty }^0\frac{\theta ({\xi }^{\prime })}{{\xi }^{\prime }-\xi }d{\xi }^{\prime }=\log q_s\text{where}{q}_s(\xi )\equiv {\left(\frac{\xi +b}{\xi +1}\right)}^{1/2}.$$

$$\mathcal{H}[\theta ](\xi )\equiv -\frac{1}{\pi }\text{vbox}-{\int }_0^{\infty }\frac{\theta ({\xi }^{\prime })}{{\xi }^{\prime }-\xi }d{\xi }^{\prime },$$

$$\log q(\xi )=\log q_s(\xi )+\mathcal{H}[\theta ](\xi ),$$

$${\mathrm{F}}^2\left[q^2\frac{dq}{d\varphi }\right]-\mathrm{T}\left[q^2\frac{d{}^2\theta }{d{\varphi }^2}+q\frac{dq}{d\varphi }\frac{d\theta }{d\varphi }\right]=-\sin \theta .$$

$$\mathcal{H}[\theta ](\varphi )=-\frac{1}{\pi }\left({\int }_{-\infty }^{-\Phi }+\text{vbox}-{\int }_{-\Phi }^{\Phi }+{\int }_{\Phi }^{\infty }\right)\frac{\theta ({\mathrm{e}}^{-{\varphi }^{\prime }})}{{\mathrm{e}}^{-{\varphi }^{\prime }}-{\mathrm{e}}^{-\varphi }}(-{\mathrm{e}}^{-{\varphi }^{\prime }})\mathrm{d}{\varphi }^{\prime },$$

$$\beta ϵq^2{q}^{\prime }-\beta \tau {ϵ}^2\left(q^2{\theta }^{\prime \prime }+q{q}^{\prime }{\theta }^{\prime }\right)+\sin \theta$$

$$\log q-\mathrm{i}\theta +\hat{\mathcal{H}}[\theta ]-\log q_s$$

$$\hat{\mathcal{H}}[\theta ](\zeta )=\frac{1}{\pi }{\int }_0^{\infty }\frac{\theta ({\xi }^{\prime })}{{\xi }^{\prime }-\zeta }\mathrm{d}{\xi }^{\prime }.$$

$$q=q_r+\bar{q}=\sum _{n=0}^{N-1}{ϵ}^n{q}_n+\bar{q}.$$

$$\bar{q}=\mathrm{i}{q}_s\bar{\theta },$$

$$\bar{q}\sim \mathcal{A}\mathcal{F}(w){\mathrm{e}}^{-\chi (w)/ϵ}.$$

$${\theta }_0$$

Full text

††volume: YYY

1 Introduction

The subject of this paper is the accurate calculation of steady-state free-surface flows where a disturbance has caused the emergence of distinct groups of gravity and capillary waves. In particular, we will discuss flow regimes that exist in the low-speed limit and numerically confirm the existence of a new class of waves predicted by Trinh & Chapman (2013b).

The canonical application of free-surface flow is the fishing-line problem first reported by Russell (1845). Figure 1 corresponds to Russell’s original illustration of a two-dimensional cross-section of the free-surface for uniform flow past a fishing line. Inspired by Russell’s report, Rayleigh (1883) modelled the fishing line as a point pressure force and demonstrated that waves of constant amplitude exist when the flow is above a certain speed (approximately $23\text{ cm/s}$ in water of infinite depth). The disturbances upstream of the fishing line are capillary waves, principally governed by surface tension, while those downstream are due to gravity (see also § 271 of Lamb 1932).

Then, in the latter half of the 20th century, it became possible to compute numerical solutions of the full nonlinear potential-flow equations for steady-state problems involving different geometries. As it pertains to the gravity-only problem, we mention as examples the works of Vanden-Broeck & Tuck (1977); King & Bloor (1987, 1990); Forbes & Schwartz (1982), although interest continues well into recent years in e.g. Părău & Vanden-Broeck (2002); Binder & Vanden-Broeck (2007); Binder et al. (2013). For reasons that will become clear, the numerical study of the types of steady-state gravity-capillary profiles illustrated in figure 1 is a much more difficult problem.

Like the methodology developed in Rayleigh (1883), theoretical analysis of such wave-structure free-surface problems usually begins by linearising the governing equations about a small parameter, say $\delta$, related to the size of the obstruction (linear geometrical theory). In addition to $\delta$, gravity-capillary flows are typically characterized by two non-dimensional parameters: the Froude number, $\mathrm{F}$ and the inverse-Bond number, $\mathrm{T}$, defined by

[TABLE]

Here, $U$ and $L$ are the chosen velocity and length scales, $\rho$ is density and $\sigma$ is the surface tension parameter. As is convention in water-wave studies, we refer to $\mathrm{T}$ directly as the Bond number. In this work, we focus on subcritical flows, i.e. those with $\mathrm{F}<1$.

As explained by Forbes (1983), linear geometrical theory indicates that for $\delta\ll 1$ there exists a critical curve, $\mathrm{T}=\mathrm{T}_{G}(\mathrm{F})$, that divides the subcritical part of the $(\mathrm{F},\mathrm{T})$ plane into two regions. Solutions with $\mathrm{T}<\mathrm{T}_{G}(\mathrm{F})$ are called Type I and are associated with the constant-amplitude waves of figure 1, while solutions with $\mathrm{T}>\mathrm{T}_{G}(\mathrm{F})$ are called Type II, and consist of localized solitary waves.

Recently, it was predicted in Trinh & Chapman (2013a, b) that several new classes of solution can occur for nonlinear geometries with $\delta=\mathcal{O}(1)$, in the limit of small Froude and Bond numbers, $\mathrm{F},\,\mathrm{T}\to 0$. The authors considered the case of a right-angled step and used techniques from exponential asymptotics to show that the typical Rayleigh (1883) bifurcation curve predicted for linear geometry widens into a band as the height of the step increases, and that within this band a range of new solutions can be found. This work aims to numerically confirm these new, Type III waves, and thus we seek to develop methods that allow us to accurately solve the governing equations for nonlinear geometries. In doing so we make use of ideas developed in Trinh (2016, 2017), where computationally problematic integral terms [cf. (20)] are removed from the governing equations through a systematic asymptotic reduction in the low-speed, $\mathrm{F},\,\mathrm{T}\to 0$ limit.

1.1 The radiation problem

Numerical solutions of steady-state free-surface problems require the use of effective radiation conditions —that no energy may come from infinity— on the edges of a truncated computational domain. However it is often unclear exactly how one can impose such a condition in practice, particularly when capillary waves are present. Stoker (1957) had previously provided several remarks on this unexpectedly challenging aspect of determining steady-state flows. As he writes (p. 175):

…it is by no means clear a priori what conditions should be imposed at infinity in order to ensure the uniqueness of a simple harmonic solution […] the steady state problem is unnatural—in the author’s view, at least—because a hypothesis is made about the motion that holds for all time, while Newtonian mechanics is basically concerned with the prediction—in a unique way, furthermore—of the motion of a mechanical system from given initial conditions.

As Stoker notes, one should in principle formulate and solve an initial-value problem where it is often sufficient to only impose boundedness of the solution at infinity. Afterwards, the time-dependent solution may then be evolved to a steady state. In practice, however, there are a myriad of reasons why solving a steady-state formulation may be necessary—beyond simple computational efficiency. For instance, the study of solutions of the direct steady-state model may yield a wealth of information about the physical or mathematical structure of the problem; information that cannot be easily obtained through time-dependent experiments.

Apart from specialized configurations (e.g. localized solitary waves), there has not yet been a proposed solution of the radiation problem, particularly for the case of Type I solutions. One approach, suggested by Grandison & Vanden-Broeck (2006) approximates the solution outside of the main computational domain using a Fourier series whose coefficients are found as part of the solution. Although this method was able to successfully compute a number of Type I solutions for flow over a small circular cylinder (non-dimensional radius $\delta=0.05$), we have found that convergence is difficult to obtain for nonlinear geometries. Other options for correctly enforcing the radiation condition include the addition of artificial viscosity, say $\mu$ [as in Părău et al. 2007] or by solving the full time-dependent problem [as in Părău et al. (2010); Moreira & Peregrine (2010)]. However, in both cases, it is not clear that taking the limit $\mu\to 0$ or $t\to\infty$ will allow recovery of the full set of solutions when $\mu=0$ and $t=\infty$.

Our main strategy forwards is to explain how, in the limit $\mathrm{F},\mathrm{T}\to 0$, the governing equations may be reduced to a single second-order differential equation, for which it is clear how the radiation condition should be applied. We solve this reduced model as a boundary-value problem, and use it to verify the new waves predicted in Trinh & Chapman (2013b). The reduced model is derived in sections 4–5, and then solved numerically in section 6. Sections 2-3 give a simple explanation for one of the new classes of waves.

2 An explanation of the new waves via the dispersion relations

Here, we give a simple explanation of one of the new classes of gravity-capillary waves that were previously described using the more sophisticated exponential asymptotics. This is done by comparing upstream and downstream dispersion relations, and is similar in spirit to arguments in Binder & Vanden-Broeck (2007).

2.1 The upstream dispersion relation

The dispersion relation that governs linear perturbations of wavenumber $\tilde{k}$ from uniform flow of speed $U$ and depth $h$ is [cf. (2.90) in Vanden-Broeck (2010)]

[TABLE]

In most typical formulations, the dispersion relation is defined in the context of the upstream quantities, thus re-scaling the expression above, we have

[TABLE]

In (3), we have used the mean upstream channel height, $L=h_{\text{up}}$, and velocity, $U=U_{\text{up}}$, for the definitions of the Froude and Bond numbers [cf. (1)]. The wavenumber has been non-dimensionalized via

[TABLE]

Thus, for a given value of $\mathrm{T}$, the dispersion relation $\mathrm{F}^{2}=G(k)$ can be plotted in the $(k,\mathrm{F}^{2})$-plane. An example profile of $G(k)$, with $T=0.15$, is shown in figure 2.

Note that if $\mathrm{T}<1/3$, the upstream dispersion relation $G(k)$ has a minimum at $(k_{G},\mathrm{F}_{G}^{2})$. The existence of the local minimum leads to two possible types of solutions within subcritical flows ($\mathrm{F}<1$):

1. Type I:

If the Froude number lies within the band $\mathrm{F}_{G}<\mathrm{F}<1$, labeled I in the figure, then there exists two distinct real wave numbers. The smaller of the two wavenumbers corresponds to gravity waves and the larger to capillary waves. The solution consists of constant-amplitude capillary waves upstream and constant-amplitude gravity waves downstream. 2. Type II:

If, however, the Froude number is less than the minimum of $G(k)$ and in regions II and III of figure 2, i.e. $0<\mathrm{F}<\mathrm{F}_{G}$, then solutions to (3) are complex-valued. Thus the wave-trains decay in the far field.

2.2 The downstream dispersion relation

In addition to the two classes described above, the analysis of Trinh & Chapman (2013b) uses techniques from exponential asymptotics to predict several new classes of waves that occur in the low-Froude, low-Bond limit for finite-bodied (nonlinear) objects. However there is a simple explanation of one of these new classes. The key idea is that if the size of the step is sufficiently large, then the velocity and length scales downstream of the step must be re-defined. Consequently, a different dispersion relation applies downstream.

We introduce the height and velocity ratios,

[TABLE]

and the downstream Froude and Bond numbers, as well as downstream wavenumbers,

[TABLE]

where by mass conservation $h_{\text{d/u}}=1/U_{\text{d/u}}$.

Applying (2) to these downstream quantities, and re-writing in terms of upstream quantities, we must have

[TABLE]

Thus in addition to the upstream dispersion relation $G(k)$, we plot the downstream dispersion relation $H(k)$ on figure 2, for the same value of $\mathrm{T}$ and with $h_{\text{d/u}}=0.8$. The downstream dispersion curve $H(k)$ has a local minimum at $(k_{H},\mathrm{F}_{H}^{2})$. Consequently, if the upstream Froude number is selected from within the band $\mathrm{F}_{H}<\mathrm{F}<\mathrm{F}_{G}$ then, although the upstream dispersion relation predicts decaying waves, the downstream dispersion relation indicates the presence of constant-amplitude waves. The result is a third class of solutions, in addition to those described in section 2.1.

1. Type III:

Decaying waves upstream and constant-amplitude gravity waves downstream. This regime is labeled as III in figure 2.

2.3 A study of the $(\mathrm{F},\mathrm{T})$-plane

The bifurcation between constant-amplitude waves and decaying waves can also be viewed in the $(\mathrm{F},\mathrm{T})$ plane, and this is shown in the top portion of figure 3. There are two important features.

Firstly, notice the thick line in the figure. This is essentially Rayleigh’s critical dispersion curve discussed in sections 1 and 2.1. That is, for $\mathrm{F}<1$ and $\mathrm{T}<1/3$, there exists a critical curve $\mathrm{T}=\mathrm{T}_{G}(\mathrm{F})$, which corresponds to the minimum over $k$ of the upstream dispersion relation $G(k;\mathrm{T})$ for each $\mathrm{T}$, where $G(k)$ is given in (3). For linear geometries, this is the only critical curve. Configurations below the curve have constant-amplitude waves (Type I solutions), and configurations above the curve have decaying waves (Type II solutions).

Secondly, notice the sequence of thinner curves in the figure. For the nonlinear step a second critical curve $\mathrm{T}=\mathrm{T}_{H}(\mathrm{F})$ exists, which corresponds to minima in the downstream dispersion relation $H(k;\mathrm{T})$, given by (7). New solutions (Type III) with decaying capillary waves and constant-amplitude gravity waves are located in the region of the $(\mathrm{F},\mathrm{T})$ plane that lies between these two curves. Note the area of this region increases with the size of the step. The thin lines in figure 3 show the downstream critical curve $\mathrm{T}_{H}(\mathrm{F})$ for various values of the depth ratio $h_{\text{d/u}}$, with the shaded area giving the region where new waves can be found for the particular case of $h_{\text{d/u}}=0.75$.

3 Asymptotics in the small Froude-Bond limit

Although our study of linear dispersion relations in the previous sections has given an insight into the new regimes in the $(\mathrm{F},\mathrm{T})$-plane, it is important to recall that (3) and (7) are derived from the study of the free-surface problem assuming a flat-bottomed channel, which is only valid sufficiently far upstream and downstream. To derive a spatially-dependent dispersion relation that connects the upstream and downstream regions, it is necessary to use the more complex exponential asymptotics, which are valid as $\mathrm{F},\mathrm{T}\to 0$. Since it is in this low-speed limit where the new waves are predicted, we re-write the dispersion relation (2) with

[TABLE]

for $\epsilon\ll 1$, following scalings from previous works (Trinh & Chapman, 2013a, b).

In the limit $\epsilon\to 0$, the wavenumber scales like $k=\hat{k}/\epsilon$. It will be shown below in section 4 that, to leading-order in $\epsilon$, the asymptotic values of the speed downstream and upstream are $\sqrt{b}$ and unity respectively, where $b$ is a parameter describing the step height. The velocity scale ratio for flow over a step is therefore

[TABLE]

In this low-speed regime, $\tanh(k)\to 1$ and the dispersion relation (3) becomes

[TABLE]

with solutions

[TABLE]

where $A=4\tau/\beta$. This is the typical small-step wavenumber in the limit of small Froude/Bond numbers, and predicts a band $A\in[0,1]$ where wavenumbers are real.

However, for a large step we must use (7) and then (6) to calculate the downstream wavenumbers as

[TABLE]

Thus, if the Froude and Bond numbers are chosen so that $A\in[1,b^{2}]$ then an additional regime is predicted, before the localised solitary waves emerge for $A>b^{2}$. The three different regions are illustrated in the lower portion of figure 3. The function $\chi(\phi)$ that gives the spatially-varying phase of the waves is found in section 3 of Trinh & Chapman (2013b). Comparison of (3.9) in that work with the expressions above shows that the wavenumbers here are just the far-field limits of $\chi$, as $|\phi|\to\infty$.

In both (11) and (12), the plus root becomes singular as $\tau\to 0$. This root therefore corresponds to capillary waves, and the minus root must correspond to gravity waves. This removes any ambiguity introduced by the extra dispersion relation, as we have exactly one upstream wavenumber that satisfies the radiation condition. The wavenumbers $k_{\text{up}}$ and $k_{\text{down}}$ can thus be used to verify that numerical solutions satisfy the radiation condition. However these wavenumbers can also be used to prescribe the behaviour of the flow at infinity, and indeed this approach will later allow us to numerically confirm the existence of the new waves in the reduced model. First, we expand on the discussion of section 1.1, and illustrate the issues surrounding the numerical solution of the full problem.

4 Boundary integral formulation

Consider steady, two-dimensional potential flow of an incompressible fluid in a channel with upstream velocity $U$, and a prescribed length scale $L$. The physical plane is shown in figure 4(a). The flow is non-dimensionalised using these characteristic scales, and we introduce complex physical coordinates, $z=x+\mathrm{i}y$, with $x$ pointing along the channel, as well as a complex potential $w=\phi+\mathrm{i}\psi$. The governing equations for the fluid are then

[TABLE]

In the kinematic condition, $\bm{n}$ denotes a normal of the channel boundary and, in Bernoulli’s equation, the curvature, $\kappa$, is defined to be positive when the centre of curvature lies inside the fluid. The derivation of this standard governing formulation is given in e.g. Chap. 2 of Vanden-Broeck (2010).

In the analysis of the potential plane, it is convenient to choose the length scale to be

[TABLE]

where $h_{\text{up}}$ is the upstream height of the channel, and thus the non-dimensional height is set to $\pi$. We also non-dimensionalise the velocity according to its upstream value $U$, so that the non-dimensional flux is $\pi$. Setting $\psi=0$ on the free surface, it must then follow that $\psi=-\pi$ on the channel bottom and consequently, the flow lies within a strip PQRS in the $w$-plane shown in figure 4(b).

It is convenient to further map the strip-flow in the $w$-plane to the upper half-$\zeta$-plane using the conformal map

[TABLE]

Thus, the free surface and channel bottom are both mapped to $\eta=0$, with $\xi>0$ the free surface and $\xi<0$ the channel bottom. The flow in the upper half-$\zeta$-plane is shown in figure 4(c), with the solid and fluid boundaries, marked PQRS, lying along the real axis.

Next, we introduce the complex velocity, $\frac{\operatorname{d\!}{}w}{\operatorname{d\!}{}z}$ by writing

[TABLE]

for fluid speed $q$ and streamline angle, $\theta$. Below, we shall recast the original system (13) in terms of $q$ and $\theta$, either written as functions of $w$ or $\zeta$.

Following the standard procedure, the requirement that $\phi$ satisfies Laplace’s equation (13a) is equivalently re-stated as a boundary-integral relationship between $q$ and $\theta$ applied on the fluid and solid surfaces. For points on the free-surface, $\xi\geq 0$, this allows us to write [cf. eqn (3.4) in Trinh & Chapman 2013a]

[TABLE]

Above, the first integral on the right hand-side is known (since we assume that the channel geometry is specified via $\theta$ for $\xi\leq 0$). The second integral is the Hilbert transform applied to the free surface, and the dash across the integral sign indicates a principal value.

In this work, we consider for convenience the case of flow over a right-angled step, given by specifying the streamline angles,

[TABLE]

where $b>1$ and $\xi=-b$ is the image of the stagnation point of the step in the $\zeta$-plane. Substitution of (18) into the first integral on the right hand-side of (17) and integrating yields

[TABLE]

In (19), we have defined the important ‘shape-function’, $q_{s}$, which encodes the channel-bottom geometry and plays a critical role in the work that follows.

Returning to the second integral on the right hand-side of (17), we shall also define the Hilbert transform operator $\mathscr{H}$, defined on the free surface $\xi\geq 0$ by

[TABLE]

so that the boundary integral equation (17) can be written compactly as

[TABLE]

To close the system, Bernoulli’s equation in (13c) is written in terms of $q$ and $\theta$ using (16). Thus, along the free surface $\psi=0$, we have

[TABLE]

The object now is to solve the integro-differential system given by the boundary integral (21a) and Bernoulli’s equation (21b) for the unknowns, $q$ and $\theta$, along the free surface given by $\xi\geq 0$ or equivalently $-\infty<\phi<\infty$ and $\psi=0$. Note however that for computational purposes we must write (20) as

[TABLE]

where $\pm\Phi$ are the limits of the computational domain. The standard approach is to discard the first and third terms in (22), which is known to introduce errors into the solution. These errors can be reduced by instead approximating the first and third terms by asymptotic solutions in the far-field, which was done for flow over a cylinder in Grandison & Vanden-Broeck (2006). Thus the numerical treatment of the Hilbert transform is key to satisfying the radiation condition.

4.1 Challenges in solving the full problem

The main challenge that confronts us in solving (21) is illustrated by figure 5(a) which shows two numerical solutions; one without surface tension ($\mathrm{T}=0$, dashed) and one with small amount of surface tension ($\mathrm{T}=2.5\times 10^{-3}$, solid). To produce figure 5, we discretise equations (21) using finite differences and solve the resulting nonlinear system using Newton’s method. This approach is used in the majority of the works referenced above (Forbes & Schwartz, 1982; Forbes, 1983; King & Bloor, 1987) and is described in detail in Vanden-Broeck (2010).

Figure 5(b) shows the Fourier spectra of the solutions $q(\phi)$, divided into either upstream or downstream components. Spectra (i) and (ii) correspond to $\mathrm{T}=0$, while spectra (iii) and (iv) correspond to $\mathrm{T}\neq 0$. In this and all subsequent figures, Fourier spectra are calculated as follows: first, the solution, $q$, is interpolated onto a regular grid and the relevant upstream ($\phi<0$) or downstream portions ($\phi>0$) extracted. The Fourier transform is taken and the results shown as a bar graph. On the (angular) wavenumber axis of each spectrum, we mark predicted gravity and capillary wavenumbers with a ‘G’ and ‘C’ as required, according to the predictions of sections 2 and 3.

There are three types of error visible in figure 5.

1. (i)

The gravity-only flow (dashed curve) should consist of constant-amplitude waves downstream and a flat (or exponentially decaying) surface upstream. However, the dotted horizontal guides show that the downstream waves are distorted over the last few periods. This distortion is primarily associated with errors in evaluating the truncated Hilbert transform (20). The transform is expressed as an integration over the whole real line, but must be truncated to a finite computational domain. 2. (ii)

Still for $\mathrm{T}=0$, the Fourier spectrum (i) shows that there are small gravity waves upstream. The radiation condition requires that only capillary waves exist upstream, so these waves are unphysical. 3. (iii)

Such errors become even more problematic for the flow with surface tension, $\mathrm{T}=2.5\times 10^{-3}$. The solid curve in figure 5(a) demonstrates the presence of spurious waves upstream, as seen by the smaller inset. For this value of $(\mathrm{F},\mathrm{T})$ the non-dimensional capillary wavenumber is $k=198$, but the upstream Fourier spectrum in (iii) shows that the dominant frequency is approximately one. Just as in the gravity-only flow, the solution does not satisfy the radiation condition. Furthermore, spectrum (iv) confirms that the downstream waves are not of the correct gravity-associated wavenumber, either.

Convergence of the numerical scheme is more difficult at larger values of $\mathrm{T}$. We emphasize that although the capillary waves shown in the smaller inset of figure 5(a) are not well resolved, the issues (i) to (iii) are not functions of the grid spacing or numerical tolerances.

Our figure 5 takes inspiration from figure 3 of Forbes (1983), who provided some initial discussion of the radiation problem for flow over a semi-circular cylinder. Forbes (1983) further indicated that without knowledge of the proper boundary conditions, it was unclear how numerical solutions could be obtained for more nonlinear flows. Similar comments on the difficulty of the radiation problem appear in the works of e.g. Forbes & Schwartz (1982), Scullen (1998), Părău & Vanden-Broeck (2002), and Grandison & Vanden-Broeck (2006) for a variety of free-surface problems.

5 Reduced models for gravity-capillary waves

The preceding discussion makes clear that one problem with numerical computation of steady-state gravity-capillary waves is the truncation of the Hilbert transform, $\mathscr{H}$, and the inability to correctly apply the radiation condition. In fact, as was argued by Trinh (2016, 2017), in the low-Froude and low-Bond limit $\epsilon\to 0$, the Hilbert transform can be systematically removed and the governing system (21) reduced. Here we present only the ideas that are most relevant to the reduction. Detailed calculations are given in the appendix. A full comparison of reduced models for gravity-only flow is given in Trinh (2017).

5.1 Removing the Hilbert transform

For the purpose of developing the reduced model, it is important to consider three main points that emerge as a result of exponential asymptotics (Chapman & Vanden-Broeck, 2006). Firstly, in the limit $\epsilon\to 0$ free-surface waves are associated with singularities in the analytic continuation of the leading-order solution $q_{s}$, given by (19). These singularities lie in the complex plane, so it is necessary to write complex versions of the governing equations. That is, we extend $\phi$ to the complex plane, and re-write $\phi\mapsto w\in\mathbb{C}$. Similarly, we re-write $\xi\mapsto\zeta\in\mathbb{C}$. Care must be taken when defining the Hilbert transform (20) in the complex plane so as to preserve the principal-value behaviour as $\zeta$ approaches the real axis. This is accomplished through a small deformation of the path of integration about the point of evaluation, which contributes a residue term.

The complex versions of the governing equations are thus

[TABLE]

where primes $(^{\prime})$ denote differentiation in $w$ and where we have defined

[TABLE]

Now in the system (23), we have moved from the free surface, where $q,\,\phi\in\mathbb{R}$, to the complex upper-half plane, where $q$ and $\phi=w\in\mathbb{C}$, and so solutions are complex-valued. To return to the free surface, where the solution is real, we must account for the contribution from the lower-half plane. This is done by adding $q$ to its complex conjugate, $q^{*}$. Understandably, these ideas of analytic continuation are subtle, but for more details, we refer readers to §5.1 and §8.4 of Trinh (2017).

Next, in the limit $\epsilon\to 0$, we expect that the solution, $q$, can be split into a base series $q_{r}$ which describes the mean flow, and a remainder term $\bar{q}$ which describes the waves [cf. Ogilvie (1968); Chapman & Vanden-Broeck (2006)]. The base series is the regular asymptotic expansion in powers of $\epsilon$ with, say, $N$ terms, and the remainder is $\mathcal{O}(\epsilon^{N})$. Thus:

[TABLE]

We also do the same for $\theta$. The inclusion of the remainder term in (25) is necessary because using a regular expansion alone predicts a flat surface at every algebraic order. In the low-speed limit, the waves are exponentially small and must be found with specialised techniques from exponential asymptotics.

Finally, a careful analysis of the governing equations using the split solution (25) reveals that computation of the mean speed $q_{r}$ only involves the Hilbert transform of known terms, that is $\hat{\mathscr{H}}[\theta_{r}]$. Further, the remainder term $\bar{q}$, and hence the form of the waves, does not depend on the Hilbert transform of $\bar{\theta}$ at leading order. Therefore the problematic term $\hat{\mathscr{H}}[\bar{\theta}]$ can be ignored in an asymptotically consistent way in the limit as $\epsilon\to 0$. Equation (23b) reduces to the simple relationship

[TABLE]

which can be substituted into (23a) to give an ordinary-differential equation (ODE) for $\bar{q}$, which is known as a reduced model.

5.2 Choosing the truncation value $N$

The derivation of the reduced model for an arbitrary truncation value $N$ is given in the appendix, but here we focus on the case $N=2$. This particular value of $N$ is chosen because it captures the functional form of the waves (Trinh, 2017). In the limit $\epsilon\to 0$, we expect the form of the waves to be

[TABLE]

If the base series (25) is truncated at $N=1$, then the resulting $N=1$ reduced model produces a solution with the correct phase function $\chi(w)$. The $N=2$ reduced model produces a solution where both $\chi$ and $\mathcal{F}$ are correct. To obtain the constant pre-factor $\mathcal{A}$, one must truncate the base series at the optimal point $\mathcal{N}(\epsilon)$, where it is shown in Chapman & Vanden-Broeck (2006) that $\mathcal{N}\to\infty$ as $\epsilon\to 0$. The choice of $N=2$ therefore correctly predicts the form of the waves, up to a multiplicative constant $\mathcal{A}$ (Trinh, 2017).

5.3 The $N=2$ reduced model for the low-Froude low-Bond limit

The first two terms in the truncated base series (25) are:

[TABLE]

where $\hat{\mathscr{H}}[\theta_{1}]$ may be evaluated explicitly, so that all terms in (28) may be written in closed form. Using the relationship (26), we may re-write (23a) as an equation for $\bar{q}$, the $N=2$ reduced model. This last substitution involves a fair amount of algebra, and so the derivation of this equation is outlined in the appendix. The $N=2$ reduced model is given by

[TABLE]

In deriving (29) we have chosen to include terms only up to $\mathcal{O}(\epsilon)$, bearing in mind the ansatz (27) which means that derivatives of $\bar{q}$ contribute a factor of $1/\epsilon$. The leading-order forcing term is $\mathcal{O}(\epsilon^{2})$. The reduced model (29) is a linear ODE with known coefficients, to be solved for $\bar{q}$. Thus (29) is computationally much more convenient than the full system (21), as the latter includes a non-local term via the Hilbert transform. Moreover, it is much simpler to apply the radiation condition to the reduced model than to the full system, and so we can be sure of selecting the unique physically-relevant solution.

Figure 6 compares numerical solutions to the full system (21) (dashed curve) and the reduced model (29) (solid curve) for gravity-only flow (i.e. with $\tau=0$). In this case, the radiation condition requires that the upstream surface is flat and it is sufficient to solve (29) by imposing the condition $\bar{q}=0$ at the first mesh point (only one condition is needed as the gravity-only problem is first-order). Note that the reduced model produces constant-amplitude waves downstream, whereas the boundary-integral formation does not [cf. figure 5(a)]. The reduced model therefore produces constant-amplitude waves of the correct frequency, but nevertheless does not agree with the boundary-integral solution downstream. This is due to truncation at $N=2$ terms.

5.4 Choosing the boundary conditions

For gravity-capillary flows, it is insufficient to treat (29) as an initial-value problem (i.e. ‘shooting’ from upstream). In particular, Type I solutions have constant-amplitude capillary waves that extend upstream to $-\infty$, so the values of $\bar{q}$ and $\bar{q}^{\prime}$ at the first mesh point are not known a priori. Even in regimes where the upstream waves decay, however, we found that enforcing a flat surface at the first mesh point leads to inaccurate and unphysical results. Instead, we treat (29) as a boundary-value problem with Sommerfeld boundary conditions at either end, using wavenumbers that satisfy the radiation condition. These wavenumbers were derived in section 3. Thus we solve (29) with the boundary conditions

[TABLE]

at the first and last mesh points, where $k$ is the relevant wavenumber given by (11) and (12) upstream and downstream respectively. By imposing (30) we have restricted ourselves to flows where the far-field solution is a train of linear waves (either constant-amplitude or decaying). However this is clearly not the only possibility, and indeed (30) only allows us to recover one of the four new classes of waves described in Trinh & Chapman (2013b).

6 Results

Here we present numerical results for the $N=2$ reduced model, and explore the low-Froude, low-Bond limit of the bifurcation diagram in figure 3 to validate the existence of the new (Type III) solutions described in Trinh & Chapman (2013b). With the boundary conditions (30) applied at either end of the computational domain, we are able to obtain good solutions to the $N=2$ model for a wide range of parameters. By keeping the values of $\beta$, $b$ and $\epsilon$ fixed (at 1, 2 and 0.5 respectively) and increasing $\tau$ we move up the vertical axis on figure 3 and vary the value of $A=4\tau/\beta$.

Figure 7 shows the two classical solutions that exist when the geometry is linear, but computed with $b=2$ (i.e. in a nonlinear geometry). Figure 7(a) shows a type I solution with $\tau=0.24$ and thus $A=0.96<1$, with constant-amplitude waves upstream and downstream. We see that the downstream and, via the inset, upstream solutions are both correctly resolved. The Fourier spectra for this solution are given in figure 9(a, b), and confirm that the dominant wavenumbers agree with the predictions from section 3, denoted G and C, and thus that the radiation condition is satisfied. However, the downstream solution is not completely ‘clean’ and small disturbances are visible on top of the main gravity waves, visible as a small peak in the spectrum at $k\approx 12$; this is shown with an arrow in figure 9(b) and is discussed further in section 6.1. Figure 7(b) shows a Type II solution, with $\tau=1.5$ and thus $A=6>b^{2}$. This is the solitary-wave solution from Rayleigh’s original classification, and has waves that decay in the far-field both upstream and downstream.

Figure 8 shows a Type III solution with $\tau=0.255$ and thus $A=1.02\in[1,b^{2}]$. This is one of the new waves from Trinh & Chapman (2013b), predicted to exist when the step height, $b-1$, is $\mathcal{O}(1)$ in size. There are constant-amplitude waves downstream, but upstream the waves decay away from the step. This is shown in greater detail in the inset. The region over which the waves decay is predicted to increase as $A\to 1$, and tests with various values of $\tau$ (not shown) confirm this. The downstream spectrum is shown in figure 9(c); and shows that the dominant frequency of the downstream solution is correct, again with a small secondary peak near the label C∗.

In summary, the careful application of the Sommerfeld boundary condition (30) to the $N=2$ reduced model (29) confirms the existence of one new class of waves predicted in Trinh & Chapman (2013b) for flow over a rectangular step. These solutions satisfy the radiation condition via the Sommerfeld boundary condition, which uses wavenumbers derived in section 3. This is confirmed by the Fourier spectra in figure 9. However, when $\epsilon=\mathcal{O}(1)$ these wavenumbers are no longer valid, and this leads to interference in the solutions. We will now discuss this in more detail.

6.1 The appearance of beating for larger values of $\epsilon$

Both the dispersion relation (10) and the $N=2$ reduced model (29) are derived from the full water-wave problem in the limit $\epsilon\to 0$. Thus for $\epsilon$ sufficiently small, solutions to (29) computed using (10) satisfy the radiation condition. For larger values of $\epsilon$, spurious secondary waves appear in the $N=2$ model, creating the ‘beating’ effect seen downstream in figure 7(a). The secondary waves arise due to interference between the two linearly independent solutions of the second-order ODE (29), which diverges from solutions of the full problem when $\epsilon=\mathcal{O}(1)$. The beating is an indication that the reduced model is being used outside its range of validity, as the low-speed dispersion relation (10) no longer predicts the far-field behaviour of the solution.

To confirm this, we can derive the far-field behaviour of (29) directly from the equation itself, and see that it differs from (10) when $\epsilon=\mathcal{O}(1)$. If we let $|\phi|\to\infty$, both $q_{0}=q_{s}$ and $q_{1}$ tend to constant values:

[TABLE]

so that far away from the step, (29) is a linear ODE with constant coefficients. The two linearly independent solutions have wavenumbers given by

[TABLE]

Comparing (32a) with (10), we see that $\mathcal{K}_{\pm}$ agrees with the upstream wavenumber $k_{\text{up}}$ for all $\epsilon$, but only agrees with the downstream wavenumber $k_{\text{down}}$ in the limit $\epsilon\to 0$. Thus if $k_{\text{up}}$ and $k_{\text{down}}$ are used in the Sommerfeld boundary conditions when $\epsilon=\mathcal{O}(1)$, beating is observed due to the discrepancy between $\mathcal{K}_{\pm}$ and $k_{\text{down}}$. The downstream spectra shown in figure 9(b)-(c) confirm that the secondary waves visible in figures 7(a) and 8 have the frequency predicted by (32a) (marked C* on the axis). Thus the interference visible in the solution is due to the use of $k_{\text{down}}$ in the downstream boundary condition in a regime where the far-field solution has frequency $\mathcal{K}_{+}$.

Thus the appearance of beating is a sign that the $N=2$ equation is reaching its limit of validity. Recall that both the dispersion relation (10) and the $N=2$ reduced model are valid asymptotic limits of the full problem when $\epsilon\to 0$; however, the wavenumbers $\mathcal{K}_{\pm}$ given in (32a) are valid for the $N=2$ equation at any $\epsilon$. Thus if the reduced model is to remain relevant in the context of the full problem, it is important that $\epsilon$ is kept small enough that (10) is still valid and beating does not occur.

7 Conclusions

In order to compute state-state gravity-capillary flow past nonlinear geometries, specialised methods are needed to select the solution that satisfies the radiation condition. In section 5, we showed that the full problem can be reduced to a linear ODE using techniques from exponential asymptotics to justify removing the problematic Hilbert transform term. This reduction is valid in the low-Froude, low-Bond limit, but places no restriction on the size of the step. This ODE was still found to be very sensitive in the sense that precise boundary conditions were needed to select solutions that satisfy the radiation condition, but this is computationally simpler to implement in the reduced model than in the boundary-integral framework of the full problem. The relevant boundary conditions were derived in section 2 by considering different dispersion relations up- and downstream.

With these boundary conditions, we were able to verify numerically the existence of new waves from Trinh & Chapman (2013b) for flow over a right-angled step, and to confirm that Rayleigh’s original bifurcation curve separates into a band when the size of the obstacle is too large to be linearised about. This was done in section 6, where one of the new classes of waves (decaying upstream, constant-amplitude downstream) was found to exist in an intermediate region between the two classical solutions [cf. figure 8].

7.1 Difficulties in applications to the full problem

One of the aims of developing these reduced models is to see whether their solutions can provide a starting guess for the iterative methods needed to solve the system of equations that comprises the full boundary integral problem (23).

The authors ran several numerical experiments, using solutions from the $N=2$ equation (29) as starting guesses for the full problem. Experiments were run with a large number ($\approx{2000}$) of grid-points and moderate parameter values ($b=1.7,\epsilon=0.5$) but these failed to converge to the correct solutions. Even on occasions where Newton’s method produced a result, the upstream wavenumbers were in serious error, and the accuracy is not significantly improved compared to solutions obtained using the initial guess $\theta=0$. As discussed in section 5.2, the $N=2$ model produces solutions that differ from the full problem by a multiplicative constant. This can be seen in the gravity-only case in figure 6, which compares solutions of the $N=2$ model and the full problem. Although the reduced model captures the wavelength correctly, there is a persistent error in amplitude. It is therefore not guaranteed that the two solutions would converge under Newton’s method. We do not believe that this conclusion should be regarded as disappointing per se, but rather that it illustrates an intrinsic difficulty that has not been noted before in an equally precise fashion.

The sensitivity of the reduced model highlights both the importance and the difficulty of numerically satisfying the radiation condition, and the difficulties in dealing with the full boundary integral problem provide further evidence that the asymptotic approach developed in Trinh (2017) and applied here is a useful one for the study of gravity-capillary flows. We will conclude with a brief overview of different approaches to the radiation problem.

7.2 A general view of the radiation problem

In regards to the general radiation problem, we believe there are two possible responses. The first is to design schemes that profit from some asymptotic or analytical property of the solution in order to impose effective boundary conditions on the numerical solver. This is what we have done here; first in simplifying the governing equations to a form that is asymptotically valid at low speeds, and then application of known properties derived using exponential asymptotics. The scheme of e.g. Grandison & Vanden-Broeck (2006), where the upstream solution is matched to a truncated Fourier series with unknown coefficients is also of this spirit. Our work here has clarified to what extent such asymptotically applied radiation conditions will work and to what extent are they disrupted [cf. section 6.1].

In practice, researchers may instead choose to add additional effects (physical or phenomenological) in order to regularize the boundary conditions. For instance, one can solve the time-dependent Euler equations and numerically investigate the $t\to\infty$ limit as in Părău et al. (2010). Or one can apply small artificial damping, $\mu$, so as to destroy the capillary waves upstream as was originally done in Rayleigh (1883). These procedures may be effective in answering questions regarding the physical phenomena, but it may be extremely difficult to recover certain mathematical structures of the full steady-state problem. For instance, time-dependent formulations will only recover stable configurations (unless, perhaps, time is reversed), or viscous formulations may alter the structure of solutions in an irreparable way. Simply said, it is unclear to what extent the limits $t\to\infty$ or $\mu\to 0$ are distinguished, and whether they account for singular effects in the entire space of steady-state solutions.

Appendix A Deriving the $N=2$ reduced model

Here, we derive the $N=2$ reduced model given by (29). We first outline the argument from Trinh (2017), which justifies the removal of the Hilbert transform for the case of gravity-only flow. Then we will describe the algebra that leads to the reduced model for gravity-capillary flow, and include a Mathematica routine that allows for the derivation of higher-order reduced models.

A.1 Removing the Hilbert transform: outline

The analysis required to justify the removal of the Hilbert transform can be easier performed on the gravity-only problem where the equations are first order. Thus we set $\tau=0$ in the complex governing equations (23), and take $\beta=1$ without loss of generality. The dependent variables $q(w)$ and $\theta(w)$ are then split into an $N$-term base series and a remainder as in (25). The governing equations are then

[TABLE]

These may be combined to give the integro-differential equation

[TABLE]

Expanding coefficients in powers of $\epsilon$, and using the fact that $\theta_{0}=0$ and $q_{0}=q_{s}$, equation (35) becomes [cf. (7.2a) in Trinh (2017)]

[TABLE]

The forcing term $R$ may be written

[TABLE]

where $\mathcal{E}_{\text{bern}}$ and $\mathcal{E}_{\text{int}}$ are the remainders after the complex governing equations (23) have been satisfied to $N$ terms by the base series. Note that this corrects equation (7.2) in Trinh (2017), where $\mathcal{E}_{\text{bern}}$ was multiplied by $\cos{\theta_{r}}/q_{r}^{2}$.

Equation (36) can be solved using the method of integrating factors, giving

[TABLE]

where

[TABLE]

Here $w_{0}$ is a singularity for $q_{s}$ in the complex plane and $w^{*}$ is an arbitrary starting point for the integration.

As written, (38) is not a closed expression for $\bar{q}$ since it involves taking the Hilbert transform of the unknown $\bar{\theta}$ in $I(w)$. Note, however, that the exponent $\chi$, which governs the phase of the waves, does not depend on the Hilbert transform at all. The same is true for gravity-capillary flow. As we outline below, $\hat{\mathscr{H}}$ is asymptotically sub-dominant in (39c) and so the problem can be vastly simplified.

The function (39c) can be analysed using the method of steepest descent (see Trinh (2016) for more details) in the limit as $\epsilon\to 0$, where the integral is dominated by the contribution from the endpoints and the singularity $w_{0}$. Expansion of the integrand at the endpoints shows that the contributions here only serve to provide additional algebraic terms to the expansion and, as is much easier seen through a normal asymptotic expansion, these depend on taking the Hilbert transform of lower-order terms. Thus $\hat{\mathscr{H}}[\theta_{r}]$ is crucial to determining the later terms of the base series $q_{r}$. However, when the contribution from the singularity is analysed, we note that

[TABLE]

is small compared to the other terms in (37). This is because the singularity $w_{0}$ lies off the free surface, and so is complex valued. The denominator of the integrand is thus bounded away from zero, and further the remainder $\bar{\theta}$ is expected to be $\mathcal{O}(\epsilon^{N})$ on the free surface, where it is evaluated during the integration. Thus, in the steepest descent analysis, the singularity produces the waves without depending on $\hat{\mathscr{H}}[\bar{\theta}]$ at leading order and so these terms can be justifiably removed from the equations without compromising the exponent or the structure of the waves.

A.2 Derivation of the $N=2$ model

The reduced model used for the numerical results in section 6 is obtained by truncating the base series after $N=2$ terms. Using the substitution (25), to leading-order in $\epsilon$ the complex governing equations (23) become

[TABLE]

whence $\theta_{0}=0$ and $q_{0}=q_{s}$. At the next order,

[TABLE]

which gives $\theta_{1}$ and $q_{1}$ as in (28c) and (28d).

To derive the reduced model, we then substitute

[TABLE]

into (23a) and expand in powers of $\epsilon$, with the base series known. By construction, the leading-order forcing term is $\mathcal{O}(\epsilon^{2})$. The coefficients of $\bar{q}$ and its derivatives are also expanded in powers of $\epsilon$, and the number of terms retained reflects the ansatz (27) which indicates that every derivative of $\bar{q}$ contributes a factor of $1/\epsilon$. For ease of use, we show in Table 1 how the reduced model can be derived in the coding language Mathematica.