Exponential asymptotics for elastic and elastic-gravity waves on flow past submerged obstacles

Abstract

Linearized flow past a submerged obstacle with an elastic sheet resting on the flow surface are studied in the limit that the bending length is small compared to the obstacle depth, in two and three dimensions. Gravitational effects are included in the two-dimensional geometry, but absent in the three-dimensional geometry; the Froude number is chosen so that gravitational and elastic restoring forces are comparable in size. In each of these problems, the waves are exponentially small in the asymptotic limit, and can be computed using exponential asymptotic methods. In the two-dimensional problem, flow past a submerged step is considered. It is found that the relative strengths of the gravitational and elastic restoring forces produce two distinct classes of elastic sheet behaviour. In one parameter regime, constant-amplitude elastic waves and gravity waves extend indefinitely upstream and downstream from the obstacle. In the other parameter regime, all waves decay exponentially away from the obstacle. The equivalent nonlinear two-dimensional geometry is then studied; this asymptotic analysis predicts the existence of a third intermediate regime in which waves persist indefinitely in only one direction, depending on whether the submerged step rises or falls. In the three-dimensional geometry, it is predicted that the elastic waves extend ahead of the submerged source, decaying algebraically in space. The form of these elastic waves is computed, and validated by comparison with numerical computations of the elastic sheet behaviour.

1. Introduction

Hydroelastic waves are waves that propagate due to flexural elasticity in a sheet or membrane that is resting on a fluid region. Early motivation for studying hydroelastic waves arose due to the use of elastic sheets as a convenient model for ice sheets floating on bodies of water (Squire et al. 1988; Squire 2007). There have been a number of studies that considered moving bodies exerting pressure on an ice surface, including Davys, Hosking & Sneyd (1985), Takizawa (1985), Schulkes, Hosking & Sneyd (1987), Milinazzo, Shinbrot & Evans (1995), Squire et al. (1996), Părău & Dias (2002), Părău & Vanden-Broeck (2011) and Dinvay, Kalisch & Părău (2019), as well as free and forced waves in ice sheets in Guyenne & Părău (2014, 2017). Motivated by recently proposed applications such as the use of piezoelectric membranes on the surface of a flow to harvest energy in Domino et al. (2018), a number of recent laboratory experiments have studied the behaviour of hydroelastic waves on smaller scales; see Akcabay & Young (2012) and Ono-dit-Biot et al. (2019). Many of these theoretical and experimental studies considered both gravitational and elastic restoring forces on the elastic sheet, producing wave behaviour known as flexural-gravity waves.

A number of studies of waves that form on elastic sheets resting on flow over submerged obstacles have been performed in linear regimes. Sturova (2014, 2015a,b) studied flow under finite or semi-infinite elastic plates past a submerged body in two dimensions using Green's function methods. The behaviour of a finite elastic plate was studied for flow over an obstacle in finite depth by Tkacheva (2015). A linearized perturbation expression for the behaviour of an elastic sheet above a point source in infinite depth was obtained in Savin & Savin (2013). Linearized geometries have also been studied using computational studies, such as the analysis in Shishmarev, Khabakhpasheva & Korobkin (2019), which examined the elastic sheet strain caused by flow past a submerged dipole in a three-dimensional channel using numerical Fourier methods.

Nonlinear models have been used to study solitary or periodic hydroelastic waves in both two and three dimensions. See, for example, the computational and asymptotic studies by Forbes (1986, 1988), Marchenko & Shrira (1991), Balmforth & Craster (1999), Părău & Dias (2002), Milewski, Vanden-Broeck & Wang (2011), Vanden-Broeck & Părău (2011), Guyenne & Părău (2012, 2015), Wang, Vanden-Broeck & Milewski (2013), Gao, Wang & Vanden-Broeck (2016), Gao, Vanden-Broeck & Wang (2018) and Trichtchenko et al. (2018). This is not an exhaustive list of research in this area; for a more comprehensive review, see Părău & Vanden-Broeck (2019). Many of these studies use a nonlinear model for the elastic sheet deformation to express the surface wave behaviour in terms of an integrable equation such as the nonlinear Schrödinger equation, which possesses soliton or periodic wave solutions.

The behaviour of nonlinear flexural-gravity waves on flow past submerged obstacles in two dimensions was considered by Stepanyants & Sturova (2021), who derived a nonlinear Schrödinger equation for weakly nonlinear perturbations above a submerged dipole, and by Semenov (2021), who studied fully nonlinear waves using a conformal mapping method developed in Forbes (1982), Forbes & Schwartz (1982) and King & Bloor (1987, 1989, 1990). The problem was formulated by applying a conformal map taking the flow region, with an unknown free-surface position, into a known domain. The resultant problem is expressed in terms of a boundary integral that can be solved numerically.

Similar analyses have been performed on related geometries, including waves on internal flow interfaces under an elastic sheet by Wang et al. (2014), finite-depth shear flow under an elastic sheet by Wang, Guan & Vanden-Broeck (2020), flow in a fluid separated by an internal elastic sheet by Părău (2018), and flow contained between two elastic sheets by Blyth, Părău & Vanden-Broeck (2011). When studying nonlinear geometries, care must be taken in choosing an appropriate model for the elastic sheet; Milewski & Wang (2013) investigated the effect of different elastic sheet models on nonlinear surface waves, and demonstrated that the choice of elastic sheet model can have a significant impact on the observed behaviour in nonlinear problems.

Elastic waves have also been the subject of experimental studies by Domino et al. (2018), Ono-dit-Biot et al. (2019) and Akcabay & Young (2012). Notably, Ono-dit-Biot et al. (2019) demonstrated that elastic waves in three dimensions exhibit qualitative behaviour similar to that of capillary waves, such as those computed in Lustri, Pethiyagoda & Chapman (2019). The scaling regimes considered in this paper are comparable to laboratory set-ups such as that in Ono-dit-Biot et al. (2019), which considers the waves that form on flexible elastic membranes suspended over an inviscid fluid. Motivated by this similarity, this paper aims to apply the exponential asymptotic techniques used in two-dimensional gravity-capillary waves in Trinh & Chapman (2013a,b), and in three-dimensional capillary waves in Lustri et al. (2019), to calculate the behaviour of hydroelastic waves.

1.1. Paper outline

We first study flexural-gravity waves in a linearized geometry generated by flow over a submerged step with small height, with elastic and gravity effects scaled by a small parameter, governing the rigidity of the elastic sheet and the Froude number (or the ratio between inertial and gravitational effects). This study is motivated by previous analyses of gravity-capillary waves by Trinh & Chapman (2013a,b), which showed that the interaction between gravitational and capillary restoring forces produces a rich variety of wave behaviour compared to capillary waves in isolation. In Chapman & Vanden-Broeck (2002), which studied capillary waves in two dimensions in the small surface tension limit, it was found using exponential asymptotic methods that capillary waves propagate with constant amplitude away from the disturbance. Similar methods were used to study gravity-capillary waves in linear (Trinh & Chapman 2013a) and nonlinear (Trinh & Chapman 2013b) geometries. Even in linear geometries, these studies demonstrated that there exist interactions between gravitational and capillary effects that affect the surface wave behaviour. Subsequent studies on elastic waves in the absence of gravity by Lustri, Koens & Pethiyagoda (2020) found that elastic wave behaviour in the absence of gravity is similar to the behaviour of the capillary waves studied in Chapman & Vanden-Broeck (2002).

This formulation will allow for direct comparison with the methods of Trinh & Chapman (2013a), which demonstrated interaction effects between gravity and capillary waves in a similar scaling limit (small surface tension and small Froude number); they found that the wave behaviour changed depending on the parameter regime. In one regime, downstream gravity waves and upstream capillary waves propagate indefinitely without decay. In the other regime, the two waves decay exponentially in space away from the submerged obstacle. The purpose of this study is to determine whether a similar variety of wave behaviour is obtained by the inclusion of gravity effects into the elastic sheet geometry studied in Lustri et al. (2020) – we will identify a similar bifurcation structure in flexural-gravity waves. Finally, we will discuss the challenges required to extend this analysis to the nonlinear case, as in Trinh & Chapman (2013b), and obtain some preliminary asymptotic results.

In the second part of this study, we investigate the behaviour of hydroelastic waves in three dimensions, and compare these results with the three-dimensional capillary wave analysis in Lustri et al. (2019). We study waves on the surface of linearized flow past a submerged point source, in the limit that the elastic rigidity is small. This geometry requires a more complicated asymptotic analysis than the two-dimensional geometry, and we therefore consider only a regime in which gravitational effects may be neglected. This analysis provides a first step to a full three-dimensional flexural-gravity wave analysis; we will discuss the challenges involved in such an analysis in the conclusion of this paper. We will determine that the hydroelastic waves do possess important similarities with capillary waves from Lustri et al. (2019). The waves are absent immediately behind the obstacle, but appear as special curves on the free surface known as ‘Stokes curves’ are crossed into the upstream region.

In both of these problems, the surface waves are exponentially small in the asymptotic parameter given by the ratio between the elastic bending length and the obstacle depth. This makes the waves impossible to compute using classical asymptotic power series methods. Instead, we use exponential asymptotic methods to obtain the surface wave behaviour. Important early examples of exponential asymptotics being used to study free-surface flow over submerged obstacles are found in Chapman & Vanden-Broeck (2002, 2006), which study exponentially small capillary and gravity waves, respectively, in regimes that make use of the full nonlinear dynamic boundary condition. These results were extended to linearized and nonlinear gravity-capillary waves in two dimensions by Trinh & Chapman (2013a,b), as well as gravity and capillary waves in linearized three-dimensional geometries by Lustri & Chapman (2013, 2014) and Lustri et al. (2019). Recently, flexural waves through an elastic sheet in the absence of gravitational effects were studied in Lustri et al. (2020), using the full nonlinear boundary condition. The present study extends directly on this body of work, exploring elastic wave effects in more detail.

The layout of this paper is as follows. We begin by introducing several models that are used in existing literature to describe the behaviour of elastic sheets, and show that these models are consistent in the linearized limit considered in the present study. We then introduce briefly the exponential asymptotic method that will be used to study flexural waves generated in elastic sheets. The remainder of the paper is divided into two parts. In the first part, we calculate the behaviour of flexural-gravity waves in a linearized two-dimensional geometry, and extend our analysis to make predictions about more complicated nonlinear geometries. In the second part, we calculate the behaviour of hydroelastic waves in a linearized three-dimensional geometry. The paper ends with conclusions and a discussion of the results, including an outline of the challenges expected in extending these results to nonlinear geometries, or introducing gravity into the three-dimensional geometry.

1.2. Waves on elastic sheets

1.2.1. Two-dimensional geometries

Elastic sheets in two dimensions have often been studied using the Cosserat model, such that the pressure jump across the sheet is related to the curvature through

(1.1)\begin{equation} p = D \left(\kappa_{ss} + \tfrac{1}{2}\kappa^3\right), \end{equation}

where $p$ is the pressure, $s$ is an arc length parameter, and $\kappa$ is the signed curvature, positive if the centre of curvature lies in the fluid region. The flexural rigidity coefficient $D$ is given by $D = Eh^3/(12(1-\nu ^2))$, where $E$ is the Young's modulus, $h$ is the plate thickness, and $\nu$ is the Poisson ratio.

By linearizing two-dimensional flows and searching for waves of the form $\exp ({\mathrm {i}(k x- \omega t)})$, it is possible to calculate the dispersion relation for flexural-gravity waves in the linear limit for flow on finite depth $L$; see, for example, the discussion in Gao et al. (2018). The dispersion relation is given by

(1.2)\begin{equation} \omega^2 = \left(g k + \frac{D k^5}{\rho}\right)\tanh(L k), \end{equation}

where $\omega$ is the angular frequency, $k$ is the wavenumber, $D$ is the flexural rigidity, and $\rho$ is the fluid density. If the depth is taken to be large, so that $\tanh (L k) \to 1$, then it is possible to determine a useful length scale for flexural-gravity waves. The phase velocity $c = \omega /k$ in this regime is given by

(1.3)\begin{equation} c^2 = \frac{g}{k} + \frac{D k^3}{\rho}. \end{equation}

The phase velocity has a minimum value $c_{{min}}$ at $k = k_{{crit}}$, where the group velocity and phase velocity are equal. These values are given by

(1.4a,b)\begin{equation} k_{{crit}} = \left(\frac{\rho g}{3D}\right)^{1/4},\quad c_{{min}} = 4 \left(\frac{Dg^3}{3\rho}\right)^{1/4}. \end{equation}

If the waves are on a steady flow past an obstacle, then waves can form only if the flow velocity $U$ exceeds $c_{{min}}$. If $U > c_{{min}}$, then the dispersion relation (1.3) gives two solutions for the wavenumber. The larger solution ($k > k_{{crit}}$) describes downstream gravity waves, and the smaller solution ($k < k_{{crit}}$) describes upstream flexural waves. From the form of $k_{crit}$ in (1.4a,b), we see that the characteristic length scale at which the waves transition from the elastic regime to the gravity regime is proportional to $l_D = (D/(\rho g))^{1/4}$, where $l_D$ is known as the ‘bending length’.

The elastic sheets used in Ono-dit-Biot et al. (2019) have values of $D$ in the range $D \approx 6 \times 10^{-6}$–$2 \times 10^{-8}\ {\rm N}\ {\rm m}^{-2}$. If such an elastic sheet is suspended above water, such that $g \approx 10\ {\rm m}\ {\rm s}^{-2}$ and $\rho \approx 10^3\ {\rm kg}\ {\rm m}^{-3}$, then the elastic bending length lies in the range $l_D \approx 0.5\times 10^{-3}$–$1\times 10^{-3}$ m. In DiMarco et al. (1993), ice sheets were calculated to have values of $D$ lying in the range $D \approx 6\times 10^9$–$9\times 10^9\ {\rm N}\ {\rm m}^{-2}$. Using the same approximate values for $g$ and $\rho$, we find that the bending length of these ice sheets is $l_D \approx 28$–30 m.

Section 2 of this study first considers linearized waves in a finite-depth channel containing a step with upstream depth $L$. The step height in the mapped potential plane after non-dimensionalization by $L$, denoted $w = \phi + \mathrm {i} \psi$, is given by $\delta$. We assume that $0 < \delta \ll 1$, producing a linearized regime. This assumption is equivalent to linearizing around small step height, or setting the ratio between the step height and channel depth – which is ${O}(\delta )$ as $\delta \to 0$ – to be asymptotically small.

After linearizing about the small step height, we then introduce a second small parameter into the problem, which describes the bending length to channel depth ratio $l_D/L$. Using the physical parameters for elastic sheets defined above, if an elastic sheet is suspended above a fluid with depth 1 cm, then it has $l_D/L \approx 0.05$–$0.2$. If an ice sheet is suspended above a channel of depth 100 m, then it has $l_D/L \approx 0.3$. Motivated by examples such as these, we are interested in studying problems in the asymptotic limit that $l_D/L$ is small.

In order to capture interactions between gravitational and elastic effects, we also set the Froude number, denoted $F$, to be small, with the relative scale chosen such that gravitational and elastic restoring effects are comparable in size. We define a new small parameter $\epsilon$ such that

(1.5a,b)\begin{equation} F^2 = \frac{U}{g L} = \beta \epsilon, \quad \left(\frac{l_D}{L}\right)^4 = \frac{D}{\rho g L^4} = \beta \tau \epsilon^4, \end{equation}

where $\beta$ and $\tau$ are chosen so that we can adjust the relationship between the Froude number $F$ and the bending length to channel ratio $l_D/L$. This particular form for the small parameter $\epsilon$ is selected so that the subsequent analysis is analogous to that of Trinh & Chapman (2013a), and the scaling of each quantity relative to $\epsilon$ is chosen so that both gravitational and elastic effects are described by our asymptotic results.

Note that the linearization step occurs before $\epsilon$ is defined. This implies that $0 < \delta \ll \epsilon \ll 1$ in the linearized problem, as we consider a full expansion of $\epsilon$ but only the leading-order equations for $\delta$. This regime allows us to establish the feasibility of the method. In § 2.5, we will study nonlinear flow over a step that is not an asymptotically small parameter. The problem formulation in this geometry will contain only the small parameter $\epsilon$, and the analysis will therefore be applicable to regimes where $0 < \epsilon \ll 1$. Much of the linearized analysis in § 2 is performed in a manner that generalizes to the nonlinear problem.

1.2.2. Three-dimensional geometries

A number of models have been used to describe elastic sheets resting on a fluid in three dimensions. The simplest linear elastic model for three dimensions is the biharmonic model. The pressure on the elastic sheet is derived using linearized beam theory, giving

(1.6)\begin{equation} p = D \varDelta^2 \eta,\end{equation}

where $\eta$ is the free-surface height, and $\varDelta$ is the biharmonic operator. This model has been used in Squire et al. (1996) and Squire (2007) to study the deformation of an elastic sheet resting on a fluid. Nonlinear approaches have been considered in the literature, such as a model based on the Cosserat theory of elastic shells, presented in Milewski & Wang (2013), Guyenne & Părău (2014) and Trichtchenko et al. (2018). This model is given in Milewski & Wang (2013) by

(1.7)\begin{align} p &= D\left(2\,\frac{\partial }{\partial x}[S(\eta_y \eta{xy} - \eta_x\eta_{yy})] + 2\,\frac{\partial }{\partial y}[S(\eta_x \eta{xy} - \eta_y\eta_{xx})] + \frac{\partial ^2}{\partial x^2}[S(1+\eta_y^2)] \right.\nonumber\\ &\quad - 2\,\frac{\partial ^2}{\partial x\,\partial y}[S\eta_x\eta_y] +\frac{\partial ^2}{\partial y^2}[S(1+\eta_x)^2] + \frac{5}{2}\,\frac{\partial }{\partial x}[S^2(1+|\boldsymbol{\nabla}\eta|^2)^{3/2}\eta_x] \nonumber\\ &\quad \left.{}+ \frac{5}{2}\,\frac{\partial }{\partial y}[S^2(1+|\boldsymbol{\nabla}\eta|^2)^{3/2}\eta_y] \right), \end{align}

where

(1.8)\begin{equation} S = \frac{(1+\eta_x^2)\eta_{yy} + (1 + \eta_y^2)\eta_{xx} - 2\eta_x\eta_y\eta_{xy}}{(1+|\boldsymbol{\nabla}\eta|^2)^{5/2}}. \end{equation}

In our study of hydroelastic waves in three dimensions, we will apply the scaling $\eta = \delta \tilde {\eta }$, where $0 < \delta \ll 1$ and $\delta$ measures the strength of the submerged source. The Cosserat model in (1.7) reduces to the biharmonic model in (1.6) under this linearization. Hence we will use the biharmonic model directly in our analysis of hydroelastic waves in three dimensions, noting that it describes behaviour produced by the commonly used nonlinear Cosserat model in the linearized regime.

In the analysis of the three-dimensional problem, we will introduce a small parameter $\epsilon$ such that $\epsilon ^3 = D/(\rho U^2 L^3)$, where $U$ is the upstream flow velocity, and $L$ is a representative length scale in the problem. The asymptotic analysis is performed on the linearized three-dimensional problem in the limit that $\epsilon$ is small. The regime in which the small-$\epsilon$ analysis performed on the linearized geometry is valid is given by $0 < \delta \ll \epsilon \ll 1$.

1.3. Exponential asymptotics

In order to study the behaviour of hydroelastic waves, we will adapt the methodology of Trinh & Chapman (2013a) for the problem of two-dimensional flexural-gravity waves, and Lustri et al. (2020) for the study of purely elastic waves in a three-dimensional setting. These studies considered waves on a free surface due to gravity or capillary effects that were exponentially small in the small Froude number and surface tension limits, respectively. In the present study, we will be performing an exponential asymptotic analysis in the limit that $\epsilon \to 0$, where $\epsilon$ is defined in (1.5a,b). In this case, it governs both gravitational and elastic effects. The limit corresponds to geometries with small Froude number, as seen previously in the exponential asymptotic study of gravity waves in Chapman & Vanden-Broeck (2006), and the ratio between the elastic bending length and the obstacle depth being small. These problems have the common property that they are singularly perturbed in the small parameter $\epsilon$, which appears in front of the leading derivative terms in the Bernoulli equation, shown after rescaling in (2.5). In singularly perturbed problems such as these, oscillatory behaviour in the solution, such as surface waves, has amplitude that is exponentially small in the asymptotic limit.

Solutions of singularly perturbed differential equations containing multiple exponential terms in the complex plane typically contain curves along which the behaviour of a subdominant exponential changes rapidly. These curves were first identified in Stokes (1864) and are known as ‘Stokes curves’. This rapid change causes exponentially small oscillations to appear in the solution, such as the elastic waves considered in the present study. Asymptotic techniques have been developed for studying this exponentially small behaviour, collectively known as ‘exponential asymptotics’. A broad summary of these techniques may be found in Boyd (1999). This investigation will apply the technique developed by Olde Daalhuis et al. (1995) and extended by Chapman, King & Adams (1998), which utilizes the rapid variation near Stokes curves in order to study exponentially small behaviour in solutions of ordinary and partial differential equations in Chapman & Mortimer (2005).

The first step in this technique is to express the solution to the differential equation system as asymptotic power series in the independent variable $w$, such as

(1.9a,b)\begin{equation} q(w) \sim \sum_{n=0}^{\infty} \epsilon^n\,q^{(n)}(w) \quad \mathrm{and} \quad \theta(w) \sim \sum_{n=0}^{\infty} \epsilon^n\,\theta^{(n)}(w) \quad \mathrm{as} \ \epsilon \rightarrow 0. \end{equation}

The series is typically divergent for singularly perturbed problems, and will not describe the behaviour of the surface waves, no matter how many series terms $q^{(n)}$ and $\theta ^{(n)}$ are calculated. This is because the waves are exponentially small in the asymptotic limit, and therefore smaller than any algebraic power of $\epsilon$. Instead, we minimize the error of the divergent series approximation by truncating the series after a particular finite number of terms. To truncate the series optimally, we follow the heuristic described in Boyd (1999) and truncate the series after its smallest term. The optimal truncation point, denoted $N_{{opt}}$, typically becomes large in the asymptotic limit, hence identifying the asymptotic form of the ‘late-order terms’ of the series (that is, the form of $q^{(n)}$ and $\theta ^{(n)}$ in the limit that $n \rightarrow \infty$) is sufficient to determine the smallest term in the series, and thus truncate the series optimally (see Chapman et al. 1998).

Dingle (1973) identified that successive terms in a divergent asymptotic series expansion generated by singularly perturbed equations like (2.5) are typically obtained by repeated differentiation of earlier terms in the series. This can be seen in the series recurrence relation (2.18), in which the series terms $q^{(n-1)}$ and $\theta ^{(n-4)}$ are differentiated once and four times, respectively. This repeated differentiation will cause singularities in earlier terms to grow in strength as $n$ increases, and therefore persist into later terms. As these singularities are differentiated repeatedly, the series terms typically diverge as the ratio between a factorial and the increasing power of a function $\chi$ that is zero at the singularity, ensuring that the late-order terms are also singular at this point. Chapman et al. (1998) proposed that the terms of a divergent asymptotic series generated in this fashion have asymptotic behaviour given by the sum of factorial-over-power ansatz expressions, each associated with a different early-order singularity. For our two-dimensional problems, the proposed factorial-over-power behaviour is given by

(1.10a,b)\begin{equation} q^{(n)} \sim \frac{Q\,\varGamma(n + \gamma)}{\chi^{n+\gamma}}\quad \mathrm{and} \quad \theta^{(n)} \sim \frac{\varTheta\,\varGamma(n + \gamma)}{\chi^{n + \gamma}} \quad \mathrm{as} \ n \to \infty, \end{equation}

where $\varGamma$ is the gamma function defined in Abramowitz & Stegun (1972), $Q$, $\varTheta$, $\gamma$ and $\chi$ are functions of $w$ that do not depend on $n$, and $\chi = 0$ at singularities of early series terms. The global behaviour of the functions $Q$, $\varTheta$, $\gamma$ and $\chi$ may be found by substituting this ansatz directly into the equations governing the terms of the asymptotic series, and matching to a rescaled local expansion of the solution in the neighbourhood of the singularity.

The late-order-term behaviour given in (1.10a,b) is related to applying a WKB (or Liouville–Green) ansatz of the form $A(w)\exp ({-\chi (w)/\epsilon })$ to the equations for $q$ and $\theta$ linearized about the truncated expansion. It is clear from this expression that $\chi$, or the ‘singulant’, determines the scaling of the exponentially small terms. The rapid change in exponentially small behaviour, or ‘Stokes switching’, occurs across curves where the switching exponential is maximally subdominant compared to the leading-order behaviour (see Dingle 1973). These curves satisfy the condition that the singulant is purely real and positive, giving the following condition that may be used to determine the possible location of Stokes lines:

(1.11a,b)\begin{equation} \mathrm{Re}(\chi) > 0,\quad \mathrm{Im}(\chi) = 0. \end{equation}

Asymptotic solutions also contain important curves known as anti-Stokes lines. These are curves that divide the complex plane into regions in which a particular exponential contribution is asymptotically small, and regions in which the exponential contribution is asymptotically large. From the WKB ansatz of the exponential contribution, it can be seen that anti-Stokes curves satisfy

(1.12)\begin{equation} \mathrm{Re}(\chi) = 0. \end{equation}

Truncating the infinite series (1.9a,b) optimally after $N$ terms gives

(1.13a,b)\begin{equation} q(w) = \sum_{n=0}^{N_{{opt}}-1} \epsilon^n\,q^{(n)}(w) + R_N(w) \quad \mathrm{and} \quad \theta(w) = \sum_{n=0}^{N_{{opt}}-1} \epsilon^n\,\theta^{(n)}(w) + S_N(w) \quad \mathrm{as} \ \epsilon \rightarrow 0, \end{equation}

where $R_N$ and $S_N$ are the exponentially small remainder terms for $q$ and $\theta$ obtained after truncation. Note that these expressions are equalities, rather than asymptotic relations. Therefore, $R_N$ and $S_N$ represent the difference between the true solution and the optimally truncated series.

The final step of the method described in Olde Daalhuis et al. (1995) requires substituting the truncated series expression back into the original problem to produce an equation for the remainder term. This remainder equation is then solved in the neighbourhood of Stokes curves, which are found using the condition in (1.11a,b). Condition (1.11a,b) is not strictly required for this step, as the location of the Stokes curves can be obtained directly from the remainder equations using late-order terms. We will use this condition, as it allows for the Stokes curves to be identified once $\chi$ has been calculated, rather than later in the analysis. This analysis shows that the exponentially small remainder that switches across the Stokes line generated by the truncated divergent series (1.13a,b) generally takes the form

(1.14a,b)\begin{equation} R_N \sim \mathcal{S} Q\,\mathrm{e}^{-\chi/\epsilon} \quad \mathrm{and} \quad S_N \sim \mathcal{S} \varTheta\,\mathrm{e}^{-\chi/\epsilon} \quad \mathrm{as} \ \epsilon \rightarrow 0, \end{equation}

where $\mathcal {S}$ is a function of $w$ that is essentially constant away from the Stokes curve, but varies rapidly in the neighbourhood of the Stokes curve. This emphasizes the important role played by the singulant in determining the behaviour of the oscillations. Importantly, if $\chi '$ is purely imaginary, then these terms correspond to exponentially small oscillations as $\epsilon \to 0$ that do not decay exponentially in space. If $Q$ and $\varTheta$ do not decay, as is the case for the oscillations in the present study, then these terms produce a train of waves with constant amplitude.

2. Two-dimensional elastic-gravity waves

2.1. Formulation

We consider a two-dimensional incompressible, irrotational, inviscid flow through a channel of finite depth over a submerged step. The upstream channel depth is given by $L$, and the upstream flow velocity is given by $U$. An elastic sheet with flexural rigidity $D$ rests on the surface of the flow. The position of the elastic sheet is denoted as $\xi (x)$. A schematic of this flow behaviour is shown in figure 1. We now non-dimensionalize the lengths of the system by the upstream depth $L$, and the velocities by the upstream flow velocity $U$.

Figure 1. Mapping of a fluid domain underneath an elastic sheet to a fixed known region of the complex potential plane. The elastic sheet is shown as a dashed line, and the rigid base is shown as an unbroken line. Steady flow follows streamlines, which are curves in the potential plane with constant $\psi$. The free surface maps to the top streamline, typically labelled $\psi = 0$, while the lower boundary typically maps to $\psi = -{\rm \pi}$; consequently, the flow region is known completely. The fluid velocity is singular at the points labelled $B$ and $C$, shown as black circles. (a) Physical fluid domain, $z=x+\mathrm {i} y$. (b) Complex potential domain, $w=\phi +\mathrm {i}\psi$.

The fluid potential satisfies Laplace's equation

(2.1)\begin{equation} \nabla^2 \phi = 0. \end{equation}

As the flow is steady, we apply a kinematic boundary condition on all boundaries,

(2.2)\begin{equation} \frac{\partial \phi}{\partial n} = 0, \end{equation}

where $n$ is the unit normal direction. On the free surface, we have the dynamic boundary condition, obtained from the Bernoulli equation,

(2.3)\begin{equation} \frac{F^2}{2}\,(|\boldsymbol{\nabla} \phi^2| - 1) + y + \frac{D}{\rho g L^4} \left(\kappa_{ss} + \frac{1}{2}\,\kappa^3\right) = 0\quad \mathrm{on} \ y = \xi(x), \end{equation}

where $\kappa$ is the curvature, defined to be positive if the centre of curvature lies within the fluid, $s$ is the arc length along the surface after non-dimensionalization, $\rho$ is the density of the fluid, $g$ is the acceleration due to gravity, $F$ is the Froude number defined in (1.5a,b), and $D$ is the flexural rigidity of the plate. The upstream flow is uniform with velocity $U$, such that

(2.4)\begin{equation} (\phi_x, \phi_y) \to (1,0) \quad \mathrm{as} \ x \to -\infty. \end{equation}

Using the quantities defined in (1.5a,b), we may rewrite (2.3) as

(2.5)\begin{equation} \frac{\beta \epsilon}{2}\,(|\boldsymbol{\nabla} \phi|^2 - 1) + y + \beta \tau \epsilon^4 \left(\kappa_{ss} + \frac{1}{2}\,\kappa^3\right) = 0\quad \mathrm{on} \ y = \xi(x), \end{equation}

where $\epsilon$ is a small parameter, and $\beta$ and $\tau$ determine the ratio between the Froude number $F$ and the elastic length ratio $l_D/L$. We note that the $\tau = 0$ problem corresponds to pure gravity waves, studied in Chapman & Vanden-Broeck (2006), while $\beta \to 0$ and $\tau = 1/\beta$ corresponds to pure elastic waves, studied in Lustri et al. (2020). In § 2.3, we will find our results to be consistent with these prior studies.

Differentiating (2.5) with respect to $s$ gives

(2.6)\begin{equation} \beta\epsilon q\,\frac{\mathrm{d} q}{\mathrm{d} s} + \frac{\mathrm{d} y}{\mathrm{d} s} + \beta \tau \epsilon^4 \left(\kappa_{sss} + \frac{3}{2}\,\kappa_s \kappa^2\right) = 0\quad \mathrm{on} \ y = \xi(x). \end{equation}

We define a complex potential $w = \phi + \mathrm {i} \psi$, where $\phi$ is the fluid potential, and $\psi$ is the streamfunction. This maps the fluid region to an infinite strip bounded by $\psi = -{\rm \pi}$ and $\psi = 0$. Noting that

(2.7a–c)\begin{equation} \kappa = \frac{\mathrm{d} \theta}{\mathrm{d} s},\quad \frac{\mathrm{d} y}{\mathrm{d} s} = \sin\theta, \quad \frac{\mathrm{d} }{\mathrm{d} s} = q\,\frac{\mathrm{d} }{\mathrm{d} \phi}, \end{equation}

we write the Bernoulli condition (2.6) in terms of $\phi$:

(2.8)\begin{align} &\beta\epsilon q^2\,\frac{\mathrm{d} q}{\mathrm{d} \phi} + \sin\theta + \beta \tau \epsilon^4 \left[ q \left(\frac{\mathrm{d} q}{\mathrm{d} \phi}\right)^3 \frac{\mathrm{d} \theta}{\mathrm{d} \phi} + 4q^2\,\frac{\mathrm{d} q}{\mathrm{d} \phi}\, \frac{\mathrm{d} \theta}{\mathrm{d} \phi}\,\frac{\mathrm{d} ^2 q}{\mathrm{d} \phi^2}+ 7q \left(\frac{\mathrm{d} q}{\mathrm{d} \phi}\right)^2\frac{\mathrm{d} ^2\theta}{\mathrm{d} \phi^2}\right.\nonumber\\ &\quad \left.{}+ 4 q^3\,\frac{\mathrm{d} ^2 q}{\mathrm{d} \phi^2}\,\frac{\mathrm{d} ^2\theta}{\mathrm{d} \phi^2} + q^3\,\frac{\mathrm{d} \theta}{\mathrm{d} \phi}\,\frac{\mathrm{d} ^3\theta}{\mathrm{d} \phi^3} + q^4\,\frac{\mathrm{d} ^4\theta}{\mathrm{d} \phi^4}\!+\! \frac{3}{2}\,q^3\,\frac{\mathrm{d} q}{\mathrm{d} \phi}\left(\frac{\mathrm{d} \theta}{\mathrm{d} \phi}\right)^3 \!+\! \frac{3}{2}\,q^4\left(\frac{\mathrm{d} \theta}{\mathrm{d} \phi}\right)^2\frac{\mathrm{d} ^2\theta}{\mathrm{d} \phi}\right] \!=\! 0. \end{align}

We also define the complex velocity $\mathrm {d} w/\mathrm {d} z = u - \mathrm {i} v$, written as $q\, \mathrm {e}^{-\mathrm {i}\theta }$. In this formulation, $q$ is the flow velocity at a point, and $\theta$ is the angle that the streamlines make with the horizontal axis. Analytically continuing (2.8) allows us to replace $\phi$ with the complex potential $w$ in (2.8) to obtain the analytically continued free-surface condition

(2.9)\begin{align} &\beta\epsilon q^2\,\frac{\mathrm{d} q}{\mathrm{d} w} + \sin\theta + \beta \tau \epsilon^4 \left[ q \left(\frac{\mathrm{d} q}{\mathrm{d} w}\right)^3 \frac{\mathrm{d} \theta}{\mathrm{d} w} + 4q^2\,\frac{\mathrm{d} q}{\mathrm{d} w}\,\frac{\mathrm{d} \theta}{\mathrm{d} w}\,\frac{\mathrm{d} ^2 q}{\mathrm{d} w^2}+ 7q \left(\frac{\mathrm{d} q}{\mathrm{d} w}\right)^2\frac{\mathrm{d} ^2\theta}{\mathrm{d} w^2}\right. \nonumber\\ &\quad \left. {}+ 4 q^3\,\frac{\mathrm{d} ^2 q}{\mathrm{d} w^2}\,\frac{\mathrm{d} ^2\theta}{\mathrm{d} w^2} + q^3\,\frac{\mathrm{d} \theta}{\mathrm{d} w}\,\frac{\mathrm{d} ^3\theta}{\mathrm{d} w^3} \!+\! q^4\,\frac{\mathrm{d} ^4\theta}{\mathrm{d} w^4} \!+\! \frac{3}{2}\,q^3\frac{\mathrm{d} q}{\mathrm{d} w}\left(\frac{\mathrm{d} \theta}{\mathrm{d} w}\right)^3 \!+\! \frac{3}{2}\,q^4\left(\frac{\mathrm{d} \theta}{\mathrm{d} w}\right)^2\frac{\mathrm{d} ^2\theta}{\mathrm{d} w}\right] \!=\! 0. \end{align}

We apply a conformal map $\zeta = \mathrm {e}^{-w}$ in order to map the fluid region from a strip in the complex potential plane to the upper half $\zeta$-plane. This map is illustrated in figure 2. We also define $\zeta = \xi + \mathrm {i} \eta$, where $\xi$ and $\eta$ are real quantities. Notably, the free surface maps to $\xi > 0$ and the base of the flow region maps to $\xi < 0$. In the mapped plane, we can apply Cauchy's theorem to obtain

(2.10)\begin{equation} \log q ={-} \frac{1}{\rm \pi} {{\int\hskip -1,05em -\,}_{-\infty}^{\infty}} \frac{\theta(\xi')}{\xi' - \xi}\,\mathrm{d}\xi'. \end{equation}

Analytically continuing this expression into the upper half-plane gives

(2.11)\begin{equation} q - \mathrm{i} \theta ={-}\frac{1}{\rm \pi}\int_{-\infty}^{\infty} \frac{\theta(\xi')}{\xi' - \zeta}\,\mathrm{d}\xi'. \end{equation}

We will define the base of the flow as a step, where $\theta = {\rm \pi}/2$ for $-(1+\delta ) < \zeta < -1$, and $\theta = 0$ for $\zeta < -(1+\delta )$ and $-1 < \zeta < 0$. Hence

(2.12)\begin{equation} q - \mathrm{i} \theta = \frac{1}{2}\log\left(\frac{\zeta+b}{\zeta+1}\right)-\frac{1}{\rm \pi}\int_{0}^{\infty} \frac{\theta(\xi')}{\xi' - \zeta}\,\mathrm{d}\xi'.\end{equation}

We note that the strip in the complex potential plane may also be analytically continued into the lower half $\zeta$-plane, which will produce complex conjugate behaviour. The full behaviour of the elastic sheet can be obtained by taking the sum of both the upper and lower half-plane contributions. We will not perform the lower half-plane calculations explicitly, but will instead add the appropriate complex conjugate contribution to the results of the upper half-plane analysis.

Figure 2. This schematic illustrates the effect of the mapping $w \mapsto \zeta$ between (a) the fluid potential domain $w=\phi +\mathrm {i}\psi$, and (b) the mapped domain $\zeta =\xi +\mathrm {i}\eta$. The mapping takes the fluid region to the entire upper half mapped plane. The elastic sheet $\psi = 0$ maps to the line $\xi > 0$, and the submerged boundary $\psi = -{\rm \pi}$ maps to the line $\xi < 0$. The elastic sheet is shown as a dashed line, and the rigid base is shown as an unbroken line. The singularities map to points that will be labelled $\xi = -b$ and $\xi = -1$.

This completes the governing equations. We treat (2.9) and (2.12) as the equations governing the analytically continued free surface. We will subsequently express the Bernoulli equation in terms of the mapped variable $\zeta$, but we will hold off until after the linearization step.

2.2. Linearization

We can linearize around the free stream for small step height. We set $b = 1 + \delta$ (assuming $0 < \delta \ll \epsilon$) and set $q = 1 + \delta \hat {q}$ and $\theta = \delta \hat {\theta }$. The linearized fluid equation is now given by

(2.13)\begin{equation} \hat{q} - \mathrm{i} \hat{\theta} = \frac{1}{2(\zeta+1)}-\frac{1}{\rm \pi}\int_{0}^{\infty} \frac{\hat{\theta}(\xi')}{\xi' - \zeta}\,\mathrm{d}\xi'. \end{equation}

The linearized Bernoulli equation is given by

(2.14)\begin{equation} \beta\epsilon\,\frac{\mathrm{d} \hat{q}}{\mathrm{d} w} + \hat{\theta} + \beta\tau\epsilon^4\,\frac{\mathrm{d} ^4\hat{\theta}}{\mathrm{d} w^4} = 0. \end{equation}

We could apply the mapping $\zeta = \mathrm {e}^{-w}$ to the Bernoulli equation, but the analysis is more straightforward in the complex potential plane. We have now fixed the problem so that the boundary follows a known curve ($\zeta > 0$ in the mapped plane, $\psi = 0$ in the complex potential plane). The linearization step is not necessary; we could apply this exponential asymptotic analysis to the fully nonlinear problem. In this case, the existence of both upstream and downstream waves on the free boundary mean that it is difficult to verify the results computationally. Progress on studying these systems has been made in the context of gravity-capillary waves in Jamshidi & Trinh (2020), but this remains a challenging numerical problem. We will discuss briefly the asymptotics of nonlinear geometries in § 2.5.

2.3. Exponential asymptotics

We write the series expression in the limit that $\epsilon \to 0$,

(2.15a,b)\begin{equation} \hat{q} \sim \sum_{n=0}^{\infty}\epsilon^{n}q^{(n)},\quad \hat{\theta} \sim \sum_{n=0}^{\infty}\epsilon^{n}\theta^{(n)}. \end{equation}

Note that including both gravity and elastic waves means that the form of the late-order terms requires powers of $\epsilon$ rather than $\epsilon ^3$, unlike in Lustri et al. (2020). Note that retaining the full power series in $\epsilon$ in a system that has been linearized in $\delta$ implies that we are considering the regime $0 < \delta \ll \epsilon \ll 1$.

The leading-order behaviour of the flow on the complex free surface is found by direct substitution, giving

(2.16a,b)\begin{equation} q^{(0)} = \frac{1}{2(\zeta+1)} = \frac{1}{2(\mathrm{e}^{{-}w}+1)},\quad \theta^{(0)} = 0. \end{equation}

The leading-order behaviour is singular at $w = \pm (2M+1) {\rm \pi}\mathrm {i}$, for $M \in \mathbb {Z}$. The singularities that matter are located at $M = \pm 1$. We will concentrate on the contributions due to the singularity at $w ={\rm \pi} \mathrm {i}$, adding the corresponding contributions afterwards. At higher orders, we obtain recurrence expressions for the complexified free surface

(2.17)\begin{gather} q^{(n)} - \mathrm{i} \theta^{(n)} ={-}\frac{1}{\rm \pi}\int_{0}^{\infty} \frac{\theta^{(n)}(\xi')}{\xi' - \zeta}\,\mathrm{d}\xi', \end{gather}

(2.18)\begin{gather}\beta\,\frac{\mathrm{d} q^{(n-1)}}{\mathrm{d} w} + \theta^{(n)} + \beta \tau\,\frac{\mathrm{d} ^4\theta^{(n-4)}}{\mathrm{d} w^4}= 0. \end{gather}

We are most interested in the form of the late-order terms, so we apply the late-order ansatz as $n \to \infty$ from (1.10a,b). In the first equation from (2.17), we neglect the integral expression, as it must be exponentially subdominant to the remaining terms in the expression. This simplification was used in Chapman & Vanden-Broeck (2002, 2006), and discussed in detail for the case of gravity waves past a ship in Trinh, Chapman & Vanden-Broeck (2011). A similar justification can be made here. At leading order as $n \to \infty$, we find $Q = \mathrm {i} \varTheta$, and

(2.19)\begin{equation} 1 - \mathrm{i}\beta\,\frac{\mathrm{d} \chi}{\mathrm{d} w} + \beta \tau \left(\frac{\mathrm{d} \chi}{\mathrm{d} w}\right)^4 = 0. \end{equation}

We can determine the form of the exponentially small contributions by solving the singulant equation (2.19). We recall that $\chi = 0$ at $w = {\rm \pi}\mathrm {i}$. This expression has four solutions, of the form

(2.20)\begin{equation} \chi = k_j (w + \mathrm{i} {\rm \pi}),\end{equation}

where $k_j$ for $j = 1,\ldots, 4$ depends on $\beta$ and $\tau$, but not $\zeta$. We denote the specific singulants as $\chi _j$ for $j = 1,\ldots,4$.

Continuing to the next order, corresponding to ${O}(q^{(n-1)})$ as $n \to \infty$, gives $Q$ and $\varTheta$ constant. To denote this clearly, we write $\varTheta = \varLambda$ and $Q = \mathrm {i}\varLambda$, where $\varLambda$ is constant in $w$, although it does depend on $\beta$ and $\tau$. This term must be determined by comparing the late-order terms with the inner problem in the neighbourhood of the singularity at $\zeta = -\mathrm {i}{\rm \pi}$. We perform this inner analysis in § A.1, and find that the prefactors are given by

(2.21)\begin{equation} \varLambda_j = \frac{k_j^8\left(k_j A_5 + k_j^2 A_6 + k_j^3 A_7 -\dfrac{1}{\beta\tau}\,A_8 \right)}{ \dfrac{2}{\beta\tau} - 6 k_j^4}, \end{equation}

where $A_j = \mathrm {i}^j \beta ^{j-3}(\beta ^3 + (4-j)\tau )$ for $j = 5, \ldots, 8$. The index choice is related to the analysis in § A.1.

Knowing that $\varTheta$ and $Q$ are constants, we can determine the value of $\gamma$ that is required for the late-order terms to be consistent with the leading-order behaviour near the singularity. The singularity in the leading order has strength 1, and this will increase by 1 at each iteration. Hence the strength of the singularity in the series term $q^{(n)}$ will be $n+1$, indicating that $\gamma = 1$. Consequently, we have fully determined the late-order asymptotic series terms (1.10a,b).

Using the methods of Olde Daalhuis et al. (1995) and Chapman et al. (1998), shown in detail in Trinh & Chapman (2013a,b), we may determine the behaviour of the surface waves that correspond with each of the four solutions of (2.20). This analysis is presented in § A.2, and gives the exponentially small wave contribution from $\chi _j$, which we denote as $\theta _{{exp},j}$, as

(2.22)\begin{equation} \theta_{{exp},j} \sim 2\, \mathrm{Re}\left[\left(\frac{1}{1 - 4 \tau \mathrm{i} k_j^3}\right)\frac{2{\rm \pi} k_j \varLambda_j}{\beta\epsilon}\exp({-k_j(w+\mathrm{i}{\rm \pi})/\epsilon})\right],\end{equation}

with a similar expression for $q_{{exp},j}$. The real part is obtained by taking the sum of the upper and lower half $\zeta$-plane contributions, which are complex conjugate values. We have therefore calculated the form of the waves, and can determine the regions in which they are present by studying the Stokes phenomenon in the system.

From the value of $\chi$ in (2.20) and the form of the exponential oscillations (1.14a,b), it is apparent that non-decaying wave behaviour can exist only if $k_j$ takes a purely imaginary value. In figures 3(a,c), we illustrate the solutions for $\tau = 1$ over a range of $\beta$, while in figures 3(b,d), we illustrate the solutions for $\beta = 1$ over a range of $\tau$. For fixed $\tau$, there exists some critical $\beta$, denoted $\beta _{{c}}(\tau )$, such that there are two values of $k_j$ with no real component for $\beta > \beta _{{c}}(\tau )$, and there are no values of $k_j$ that are purely imaginary if $\beta$ is less than this critical value. Conversely, for fixed $\beta$, there exists a critical value of $\tau$, denoted $\tau _{{c}}(\beta )$, such that there are two values of $k_j$ with no real component for $\tau < \tau _{{c}}(\beta )$, and no purely imaginary values of $k_j$ if $\tau$ exceeds this critical value. Hence non-decaying wave behaviour exists in the solution only if $\beta > \beta _{{c}}(\tau )$, or equivalently, $\tau < \tau _{{c}}(\beta )$. The wave behaviour in the $\beta$–$\tau$ parameter space is illustrated in figure 4.

Figure 3. Real and imaginary components of $k_1$ (red) and $k_2$ (black), $k_3$ (blue) and $k_4$ (magenta) for (a,c) $\tau = 1$, and (b,d) $\beta = 1$. Dashed lines indicate multiple solutions $k_j$ taking identical values. In (a,c), there is a critical value of $\beta$ above which $\mathrm {Re}(k_{1,2}) = 0$. In (b,d), there is a critical value of $\tau$ below which $\mathrm {Re}(k_{1,2}) = 0$.

Figure 4. Illustration of the wave behaviour as $\beta$ and $\tau$ are varied. In the unshaded region, corresponding to $\beta < \beta _{{c}}$, all of the surface waves decay spatially away from the obstacle. In the shaded region, corresponding to $\beta > \beta _{{c}}$, the surface behaviour contains non-decaying waves in the downstream and upstream directions, driven by gravitational and elastic restoring forces, respectively.

As (2.19) is a quartic equation, the four solutions may be computed exactly. We denote the solutions according to their asymptotic behaviour in the limit that $\tau \rightarrow 0$ for fixed $\beta$, which corresponds to the gravity wave limit. The four solutions have the behaviour

(2.23a–c)\begin{equation} \left.\begin{gathered} k_1 ={-}\frac{\mathrm{i}}{\beta} + \mathcal{O}(\tau^{1/6}), \quad k_2 ={-}\frac{\mathrm{i}}{\tau^{1/3}} + \frac{\mathrm{i}}{3\beta} + \mathcal{O}(\tau^{1/6}),\\ k_{3,4} ={-}\frac{\mathrm{i} \pm\sqrt{3}}{2\tau^{1/3}} + \frac{\mathrm{i}}{3\beta}+ \mathcal{O}(\tau^{1/6}). \end{gathered}\right\} \end{equation}

For each of these values of $k_j$, we denote the corresponding singulant as $\chi _j$. The exponentially small waves given by the solution $\chi _1$ tend to known gravitational wave behaviour from Chapman & Vanden-Broeck (2006) in this limit, while the amplitude of the remaining waves tends to zero.

The three wave contributions associated with $\chi _2$, $\chi _3$ and $\chi _4$ correspond to three elastic wave contributions, and are equivalent in the limit that $\beta \to 0$ and $\tau \to 1/\beta$ to those found in Lustri et al. (2019). One of these wave contributions, corresponding to the singulant $\chi _2$, produces waves that do not decay in space away from the obstacle in this limit, while the remaining wave contributions decay exponentially in space. This is consistent with the behaviour identified in Lustri et al. (2019).

The exact form of these solutions may be determined using a computational algebra package, allowing us to determine the critical values of $\beta$ and $\tau$ exactly. This bifurcation corresponds to a branch point in the exact solutions at $27\beta _{{c}}^3 - 256 t_{{c}} = 0$, or

(2.24)\begin{equation} 3\beta_{{c}} = 2^{8/3}\tau^{1/3}_{{c}}. \end{equation}

This curve in the $\beta$–$\tau$ parameter space divides solutions with non-decaying oscillations corresponding to gravity waves (corresponding to the solution $k_1$) and elastic waves (corresponding to the solution $k_2$), and solutions in which all four of the wave contributions decay. This result corresponds to setting the flow velocity $U$ to be equal to $c_{{min}}$ in (1.4a,b). It is not surprising that we recover the critical speed from the infinite-depth problem, as the linearization step requires the obstacle depth to be small compared to the channel depth. These results are therefore consistent with the phase velocity behaviour (1.3). It is impossible in this regime to choose parameters such that only one of the upstream elastic waves or downstream gravitational waves decays, while the other has constant amplitude. The solution either contains both non-decaying surface wave contributions, or neither. We will see later that this is not necessarily true for nonlinear geometries.

When non-decaying waves are present, and therefore $k_j$ is imaginary, the surface behaviour can be computed by integrating (2.22) with respect to $x$. Noting that $\phi = x$ to leading order in $\delta$, $\varLambda _j$ is real when $k_j$ is imaginary, and recalling that the wave amplitude is scaled by $\delta$, we can obtain a wave amplitude satisfying

(2.25)\begin{equation} \mathrm{amplitude} \sim \frac{2{\rm \pi} \varLambda_j \delta}{\beta(1-4\tau \mathrm{i} k_j)}\exp({-|k_j|{\rm \pi}/\epsilon}) \quad \mathrm{as} \ \epsilon \to 0. \end{equation}

From the leading-order behaviour $q^{(0)}$ in (2.16a,b), it can be determined that the non-dimensional step height is $\delta /2$ to leading order in $\epsilon$. Hence we are able to determine a relationship between the wave amplitude and the step height in the linearized regime.

From this expression, the amplitude can be expressed in terms of the upstream flow speed and the step height, rather than $\delta$ and $\epsilon$. The step height is $\delta /2$ in the limit that $\epsilon \to 0$. The upstream flow speed can be incorporated using (1.4a,b) and (1.5a,b) to give $U/c_{{min}} = 3^{1/4}/(4\tau \epsilon ^3)$, which can be solved for $\epsilon$. This allows the amplitude expression in (2.25) to be written in terms of physical properties of the flow geometry, and used to determine quantities such as the wave energy, which scales with the square of the amplitude.

2.4. Stokes structure

The Stokes structure of the analytically continued free surface is presented in figures 5(a) and 5(b) for $\beta > \beta _{{c}}$ and $\beta < \beta _{{c}}$, respectively. Schematics representing the physical flow behaviour are shown in figures 5(c) and 5(d) for $\beta > \beta _{{c}}$ and $\beta < \beta _{{c}}$, respectively.

Figure 5. Schematics of the analytically continued free surface are presented in (a,b). Stokes curves are shown as blue lines that originate at the leading-order singularity, shown as a black circle at $w = -\mathrm {i}{\rm \pi}$. The points where the Stokes curve intersect the free surface are depicted as blue circles. If the waves are constant amplitude, then the blue circles are filled. If the waves decay in space, then the circles are unfilled. The regions in which the wave contributions are present are indicated by arrows above the schematics. The cases for $\beta > \beta _{{c}}$ and $\beta < \beta {{c}}$ are shown in (a,b), respectively. In (a), the free surface contains constant-amplitude upstream elastic waves and downstream gravity waves, which are both switched on across the Stokes curve that follows $\mathrm {Re}(w) = 0$. The surface intersects two other Stokes curves, which switch on decaying free-surface waves. In (b), the free surface contains only decaying waves in the far field. This physical behaviour is shown in (c,d) for $\beta > \beta _{{c}}$ and $\beta < \beta c$, respectively. In (c), the constant-amplitude upstream and downstream waves are depicted. The decaying waves are exponentially small compared to the constant-amplitude waves as $\epsilon \to 0$, and therefore not shown. In (d), all four wave contributions are shown, and it can be seen that they decay in space.

Figure 5(a) shows the Stokes structure on the analytically continued free surface for $\beta > \beta _{{c}}$. Here, elastic and gravity wave contributions are switched on across a Stokes curve that extends vertically from the singularity at $w = -\mathrm {i} {\rm \pi}$, corresponding to $\chi _1$ and $\chi _2$, respectively. The gravity waves extend downstream from the Stokes curve, while the elastic waves extend upstream. As $k_1$ and $k_2$ are imaginary, the waves do not decay in space, but rather persist with constant amplitude.

Two other Stokes curves are present on the analytically continued free surface, corresponding to $\chi _3$ and $\chi _4$. These contributions cause rapidly decaying waves to appear downstream and upstream from the obstacle, respectively. In this case, the direction of propagation is determined by the sign of $\mathrm {Re}(k_3)$ and $\mathrm {Re}(k_4)$. As $\mathrm {Re}(k_3) > 0$, the waves must appear on the downstream side of the Stokes curve, as they would grow exponentially in the upstream direction. Conversely, as $\mathrm {Re}(k_3) < 0$, the waves must appear only in the upstream direction, where they decay exponentially in space.

A schematic of the physical behaviour of this system is shown in figure 5(c), which depicts constant-amplitude gravity and elastic waves in the downstream and upstream directions, respectively. The decaying waves are exponentially small compared to the constant-amplitude waves as $\epsilon \to 0$ even at the point in the surface where they first appear, so they are not depicted in the schematic. In this schematic, the gravitational waves are represented with a larger amplitude than the elastic waves. From figures 3(a,c), we see that $0 > \mathrm {Im}(k_1) > \mathrm {Im}(k_2)$ in the region $\beta > \beta _{{c}}$. From the form of the exponentially small terms in (2.22), this implies that the gravity waves associated with $\chi _1$ must have a greater amplitude than the elastic waves associated with $\chi _2$.

Figure 5(b) shows the Stokes structure on the analytically continued free surface for $\beta < \beta _{{c}}$. As $\mathrm {Re}(k_j) < 0$ for $j = 1,3$, the waves associated with $\chi _1$ and $\chi _3$ must decay downstream from the corresponding Stokes curve. Conversely, as $\mathrm {Re}(k_j) > 0$ for $j = 2,4$, the associated waves must decay upstream from the corresponding Stokes curve. A schematic of this physical configuration is shown in figure 5(d). All four wave contributions decay exponentially in space, meaning that the surface far upstream and downstream from the obstacle must be flat, with no waves present.

2.5. Nonlinear problem

We note that an exponential asymptotic analysis of gravity-capillary waves in nonlinear regimes in Trinh & Chapman (2013b) revealed a complicated wave structure, including second-generation Stokes switching (see Chapman & Mortimer 2005), which was caused by interactions between gravity and capillary effects. We will outline the steps required in order to analyse flexural-gravity waves in a nonlinear regime, and determine the singulant equation for general flow over topography. We will then consider the singulant behaviour for flow over a step in a nonlinear regime. This geometry again corresponds to figure 1, although the step height is no longer small.

The formulation of this problem without linearization is largely analogous to the previous analysis, except that the analytically continued nonlinear dynamic boundary condition is given by (2.9). The remaining steps follow essentially the same format. We obtain a system of recurrence equations for the series terms that is similar to (2.17)–(2.18), which can be solved to determine the algebraic series terms for the flow behaviour. The leading-order behaviour of the flow is given by

(2.26a,b)\begin{equation} \theta^{(0)} = 0,\quad q^{(0)} = \left(\frac{\mathrm{e}^{{-}w}+b}{\mathrm{e}^{{-}w} + 1}\right)^{1/2}. \end{equation}

We now pose a late-order ansatz identical to (1.10a,b) and apply this to the recurrence relation. Matching the resultant expression in the limit that $n \rightarrow \infty$ gives a singulant equation at leading order:

(2.27a,b)\begin{equation} 1 - \beta (q^{(0)})^3\,\frac{\mathrm{d} \chi}{\mathrm{d} w} + \beta \tau (q^{(0)})^4 \left(\frac{\mathrm{d} \chi}{\mathrm{d} w}\right)^4 = 0,\quad \chi(-\mathrm{i}{\rm \pi}) = 0. \end{equation}

This is a nonlinear differential equation that depends on the leading-order flow behaviour. Even without solving this differential equation, we are able to make some observations regarding the upstream and downstream flow behaviour. We denote $q^{(0)}$ in the limit that $w \rightarrow \infty$ as $q_{{down}}$, corresponding to the velocity far downstream from the obstacle. Similarly, we denote upstream flow velocity, corresponding $q^{(0)}$ in the limit $w \to -\infty$, as $q_{{up}}$.

The nonlinear system has different critical values of $\beta$ and $\tau$ for upstream and downstream waves. The downstream critical values are given by

(2.28)\begin{equation} 3 q_{{down}}^{8/3} \beta_{{c,down}} = 2^{8/3}\tau^{1/3}_{{c,down}}, \end{equation}

while the upstream critical values are given by

(2.29)\begin{equation} 3q_{{up}}^{8/3} \beta_{{c,up}} = 2^{8/3}\tau^{1/3}_{{c,up}}, \end{equation}

where the subscripts indicate whether the critical value describes the upstream or downstream region. For the step geometry in figure 1, we have $q_{{down}} = b^{1/2}$, while $q_{{up}} = 1$. This gives

(2.30a,b)\begin{equation} \beta_{{c,up}} = \left(\frac{256}{27}\,\tau_{{c,up}}\right)^{1/3},\quad \beta_{{c,down}} = \left(\frac{256}{27b^4}\,\tau_{{c,down}}\right)^{1/3}. \end{equation}

There are three possible elastic sheet behaviours. If $\beta < \beta _{{c,down}}$, then all waves on the free surface must decay in space away from the step. If $\beta > \beta _{{c,up}}$, then the surface can contain non-decaying gravitational waves downstream from the step, and non-decaying elastic waves upstream from the step. These configurations were both possible in the linearized problem. However, if $\beta _{{c,down}} < \beta < \beta _{{c,up}}$, then any downstream gravitational waves have non-decaying amplitude, while all elastic effects must decay in space away from the step. The parameter regimes are illustrated for a step with $b = 2$ in figure 6(a).

Figure 6. Possible wave behaviour as $\beta$ and $\tau$ are varied, for (a) an upwards step with $b = 2$, and (b) a downwards step with $b=2$. In (a,b), there are three regions. In the unshaded region, all elastic and gravitational wave contributions must decay in space away from the step. In the dark-shaded region, both elastic and gravitational waves can propagate indefinitely with constant amplitude. In (a), the lightly shaded region corresponds to $\beta _{{c, down}} < \beta < \beta _{{c, up}}$, and the surface can support non-decaying waves only in the downstream region, driven by gravity. Any upstream elastic waves must decay in space. In (b), the lightly shaded region corresponds to $\beta _{{c, up}} < \beta < \beta _{{c, down}}$, and the surface can support non-decaying waves only in the upstream region, driven by elastic forces. Panels (c,d) depict the one-sided wave profiles for (a,b), respectively. In (c), the schematic shows an upwards step that produces downstream gravity waves with constant amplitude, while the upstream waves decay. In (d), the schematic shows a downwards step that produces upstream elastic waves with constant amplitude, while the downstream gravity waves decay. Neither of these two behaviours is possible in the linearized system.

It is also possible to consider a downwards step, such that $\theta _0 = 0$ for $-b < \zeta < -1$. In this case, the leading-order solution is given by

(2.31a,b)\begin{equation} \theta^{(0)} = 0,\quad q^{(0)} = \left(\frac{\mathrm{e}^{{-}w}+1}{\mathrm{e}^{{-}w} + b}\right)^{1/2}. \end{equation}

The new critical values are instead given by

(2.32a,b)\begin{equation} \beta_{{c,up}} = \left(\frac{256}{27}\,\tau_{{c,up}}\right)^{1/3},\quad \beta_{{c,down}} =\left(\frac{256b^4}{27}\,\tau_{{c,down}}\right)^{1/3}. \end{equation}

If $\beta < \beta _{{c,up}}$, then all waves on the free surface must decay in space away from the step. If $\beta > \beta _{{c,down}}$, then the flow exceeds both critical values of $\beta$, and can support non-decaying gravitational waves downstream from the step, and non-decaying elastic waves upstream from the step. If $\beta _{{c,up}} < \beta < \beta _{{c,down}}$, then any downstream gravitational waves must decay in space, while non-decaying elastic effects are possible. The parameter regimes are illustrated for a downwards step with $b = 2$ in figure 6(b).

This behaviour is consistent with the dispersion relation for finite-depth flow (1.2). We denote the upstream depth of the channel as $L_{{u}}$ and the downstream depth as $L_{{d}}$; for an upwards step, $L_{{d}} < L_{{u}}$. This means that the minimum value for the phase velocity, $c_{{min}}$, differs on either side of the step. As $L$ increases, it can be seen from the dispersion relation (1.2) that the minimum value of the phase speed for flexural-gravity waves also increases. In figure 6(a), we see that the geometry permits an intermediate region in which there are only downstream gravity waves. This corresponds to the case where the flow velocity exceeds $c_{{min}}$ in the downstream region, but is lower than the higher value of $c_{{min}}$ obtained in the upstream region. The converse is true for downstream steps, where $L_{{u}} < L_{{d}}$; in this case, the minimum phase velocity for waves is greater in the downstream region than the upstream region, leading to flow geometries described in figure 6(b) that contain only upstream waves.

This is not a full analysis of the nonlinear problem. We would need to study $\chi$ in order to identify the position of Stokes curves in the problem, and therefore determine the location at which the waves appear. This requires determining the solution to (2.27a,b), which would likely necessitate a computational study. Nonlinear problems can also contain more complicated switching behaviour, such as second-generation Stokes switching; see Body, King & Tew (2005) and Chapman & Mortimer (2005). This was found to exist in some nonlinear gravity-capillary wave regimes in Trinh & Chapman (2013b). Any conclusions reached for the nonlinear system would require validation against numerical simulations, but the presence of surface waves in both directions far from the step makes it challenging to obtain sensible boundary conditions for the flow behaviour. For a detailed description of the numerical challenges involved in studying these systems, and substantial progress in overcoming these obstacles, see the numerical analysis of two-dimensional gravity-capillary waves in Jamshidi & Trinh (2020). A full analysis of the nonlinear problem is therefore beyond the scope of the present study.

An important difference between the analysis of the linear problem and any full nonlinear analysis is that the wavelength of the flexural-gravity waves will depend on $b$, and therefore the step height. This may be seen by the inclusion of $q^{(0)}$ in the singulant equation (2.27a,b). The explicit dependence can be computed by solving (2.27a,b) using the asymptotic behaviour of $q^{(0)}$ in the limit that $w \to -\infty$ for upstream waves, and $w \to \infty$ for downstream waves. This phenomenon is predicted by the dispersion relation (1.2), which depends explicitly on the depth of the channel, and is consistent with other related exponential asymptotic studies such as Chapman & Vanden-Broeck (2002, 2006), Trinh & Chapman (2013b) and Lustri et al. (2020). This wavelength selection did not occur in the linearized problem, as the flow was linearized around an unperturbed flow of constant depth; this constant depth determines the wavelength of the flexural-gravity waves to leading order in $\delta$.

3. Three-dimensional hydroelastic waves

We consider a three-dimensional incompressible, irrotational, inviscid flow of infinite depth with a submerged point source at depth $H$ and upstream flow velocity $U$. An elastic sheet with flexural rigidity $D$ rests on the surface of the flow. In three dimensions, we non-dimensionalize velocity by the upstream flow velocity $U$, and distance by a reference length scale $L$. The position of the elastic sheet is denoted as $\xi (x,y)$. The flow therefore has a non-dimensionalized source depth $h = H/L$.

A schematic of this flow geometry is shown in figure 7. A computed elastic sheet solution for $\epsilon = 0.15$, where $\epsilon ^3 = D/(\rho U^2 L^3)$, is shown in figure 8. The waves that form in the elastic sheet persist upstream from the obstacle, in a similar fashion to capillary waves. This scaling regime neglects gravitational effects; this decision is justified by experimental work such as Ono-dit-Biot et al. (2019), in which gravitational effects are present, but the gravitational wavelength is sufficiently large that they are not apparent in the experimental results. Hence we can predict the behaviour of the hydroelastic waves in such a set-up without incorporating gravitational effects.

Figure 7. Prescribed fluid configuration for three-dimensional flow with unit non-dimensionalized velocity past a source with non-dimensionalized depth $h$. The shaded region represents the position of the free surface $\xi (x,y)$, and the cross represents the position of the source. The flow region lies below the free surface, and the mean flow is moving from left to right, with flow velocity $U$ in the unscaled problem. Elastic waves form upstream from the obstacle.

Figure 8. Computed three-dimensional linearized flow for $\epsilon = 0.15$ past a source with unit depth, satisfying the system given in (3.6)–(3.9). The surface along the line $y = 0$ is shown as a black curve, with visible ripples ahead of the source, submerged one unit under $(x,y) = (0,0)$.

3.1. Governing equation

The flow is governed by Laplace's equation in three dimensions,

(3.1)\begin{equation} \nabla^2\phi = 0,\quad -\infty < z < \xi(x,y), \end{equation}

with kinematic boundary condition

(3.2)\begin{equation} \xi_x\phi_x + \xi_y\phi_y = \phi_z,\quad z = \xi(x,y). \end{equation}

As discussed in § 1.2, we apply a biharmonic model to describe the behaviour of the elastic sheet in the linear regime. This gives

(3.3)\begin{equation} \tfrac{1}{2}\left(|\boldsymbol{\nabla}\phi|^2-1\right) +\epsilon^3(\xi_{xxxx} + 2 \xi_{xxyy} + \xi_{yyyy}) = 0,\quad z = \xi(x,y), \end{equation}

where $\epsilon ^3 = D/(\rho U^2 L^3)$. This quantity corresponds to the ratio between the Froude number and bending length ratio presented in (1.5a,b). We are concerned with the free-surface behaviour in the limit $0 < \epsilon \ll 1$, corresponding to a regime in which gravity is neglected, and inertial effects are large compared to the elastic restoring force. Since the flow is uniform in the far field, $\phi _x \rightarrow 1$. The source condition is set to

(3.4)\begin{equation} \phi \sim \frac{\delta}{4{\rm \pi}\sqrt{x^2 + y^2 + (z+h)^2}} \quad \mathrm{as} \ (x,y,z) \rightarrow (0,0,-h). \end{equation}

Finally, we can prescribe that the solution satisfies a radiation condition, with waves present directly upstream from the singularity.

We are concerned with the limit $0 < \delta \ll \epsilon$, describing a weak source. In this case, the flow disturbance due to the source effect is small, and the equations may be linearized in $\delta$ about a uniform stream while retaining the full asymptotic behaviour in the small-$\epsilon$ limit. We are therefore studying the combined asymptotic parameter regime $0 < \delta \ll \epsilon \ll 1$.

3.2. Linearization

We linearize about uniform flow by setting

(3.5a,b)\begin{equation} \phi = x + \delta\tilde{\phi},\quad \xi = \delta\tilde{\xi}, \end{equation}

to give, at leading order in $\delta$,

(3.6)\begin{gather} \nabla^2\tilde{\phi} = 0 ,\quad -\infty < z < 0, \end{gather}

(3.7)\begin{gather} \tilde{\phi}_z - \tilde{\xi}_x = 0,\quad z= 0, \end{gather}

(3.8)\begin{gather} \tilde{\phi}_x -\epsilon^3 \left(\tilde{\xi}_{xxxx} +2 \tilde{\xi}_{xxyy} + \tilde{\xi}_{yyyy}\right)=0 ,\quad z = 0, \end{gather}

where the boundary conditions are now applied on the fixed surface $z = 0$. The far-field conditions imply that $\tilde {\phi } \rightarrow 0$ as $x^2 + y^2 + z^2 \rightarrow \infty$, while near the source, the singular behaviour is given by

(3.9)\begin{equation} \tilde{\phi} \sim \frac{1}{4{\rm \pi}\sqrt{x^2 + y^2 + (z+h)^2}} \quad \mathrm{as} \ (x,y,z) \rightarrow (0,0,-h). \end{equation}

3.3. Series expression

We first expand the fluid potential and free-surface position as power series in $\epsilon$,

(3.10a,b)\begin{equation} \tilde{\phi} \sim \sum_{n=0}^{\infty}\epsilon^{3n}\phi^{(n)},\quad \tilde{\xi} \sim \sum_{n=0}^{\infty}\epsilon^{3n}\xi^{(n)}, \end{equation}

to give for $n \geq 0$

(3.11)\begin{gather} \nabla^2 \phi^{(n)} = 0,\quad -\infty < z < 0, \end{gather}

(3.12)\begin{gather}\phi^{(n)}_z - \xi^{(n)}_x = 0,\quad z= 0, \end{gather}

(3.13)\begin{gather}{\phi}_x^{(n)} -{\xi}_{xxxx}^{(n-1)} -2 {\xi}_{xxyy}^{(n-1)} - {\xi}_{yyyy}^{(n-1)} =0 , \quad z = 0, \end{gather}

with the convention that $\xi ^{(-1)} = 0$. The far-field behaviour tends to zero at all orders of $n$, and the singularity condition (3.9) is applied to the leading-order expression, giving

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

The leading-order solution is given by

(3.15)\begin{gather} \phi^{(0)} = \frac{1}{4{\rm \pi}\sqrt{x^2 + y^2 + (z+h)^2}} - \frac{1}{4{\rm \pi}\sqrt{x^2 + y^2 + (z-h)^2}}, \end{gather}

(3.16)\begin{gather}\xi^{(0)} ={-}\frac{x h}{2 {\rm \pi}(y^2 + h^2)\sqrt{x^2+y^2+h^2}} - \frac{1}{2{\rm \pi}(y^2+h^2)}, \end{gather}

where the leading-order free-surface behaviour is set to be undisturbed far behind the source.

3.4. Late-order terms

In order to optimally truncate the asymptotic series prescribed in (3.10a,b), we must determine the form of the late-order terms. To accomplish this, we make a factorial-over-power ansatz with the form

(3.17a,b)\begin{equation} \phi^{(n)} \sim \frac{\varPhi(x,y,z)\,\varGamma(3n+\gamma)}{\chi(x,y,z)^{3n+\gamma}}\quad \mathrm{and} \quad \xi^{(n)} \sim \frac{\varXi(x,y)\,\varGamma(3n+\gamma)}{\chi(x,y,0)^{3n+\gamma}} \quad \mathrm{as} \ n \rightarrow \infty, \end{equation}

where $\gamma$ is a constant. In order that (3.17a,b) is the power series developed in § 3.3, we require that the singulant, $\chi$ satisfies

(3.18)\begin{equation} \chi = 0 \quad \mathrm{on} \ x^2 + y^2 + (z\pm h)^2 = 0, \end{equation}

where the sign chosen depends upon which of the two singularities is being considered. For complex values of $x$, $y$ and $z$, this defines a four-dimensional hypersurface. Irrespective of which singularity is under consideration, this hypersurface intersects the four-dimensional complexified free surface on the two-dimensional hypersurface satisfying $x^2 + y^2 + h^2 = 0$.

3.4.1. Calculating the singulant

Applying the ansatz expressions in (3.17a,b) to the governing equation (3.11) and taking the first two orders as $n \rightarrow \infty$ gives, for $z \leq 0$,

(3.19)\begin{gather} \chi_x^2 + \chi_y^2 + \chi_z^2 = 0, \end{gather}

(3.20)\begin{gather}2\varPhi_x\chi_x + 2\varPhi_y\chi_y + 2\varPhi_z\chi_z ={-}(\chi_{xx}+\chi_{yy}+\chi_{zz}), \end{gather}

while the boundary conditions on $z = 0$ at leading order become

(3.21)\begin{gather} -\chi_z\varPhi + \chi_x \varXi = 0, \end{gather}

(3.22)\begin{gather}\chi_x\varPhi + (\chi_x^4 +2 \chi_x^2\chi_y^2 + \chi_y^4)\varXi = 0. \end{gather}

The system in (3.21)–(3.22) must have non-zero solutions, which requires

(3.23a,b)\begin{equation} \chi_z ={-}\frac{\chi_x^2}{(\chi_x^2+\chi_y^2)^2}, \quad \varXi ={-}\frac{\chi_x}{\chi_x^2+\chi_y^2}\,\varPhi. \end{equation}

Applying (3.23a,b) to (3.19) evaluated on $z=0$ gives a singulant equation for $\chi$ on the free surface:

(3.24)\begin{equation} \chi_x^4 + \left(\chi_x^2 + \chi_y^2\right)^5 = 0. \end{equation}

This expression is similar to the capillary wave singulant equation from Lustri et al. (2019), with a different power in the second term. The subsequent analysis is therefore similar, and we include an outline of the details. Because the singularity lies below the fluid surface, we must solve (3.24) for complex $x$ and $y$ with the boundary condition

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

Solving (3.24)–(3.25) using Charpit's method gives

(3.26)\begin{equation} \chi ={\pm}\frac{3 z_j h^{1/3}s^{5/3}(s-x)}{2h^2+5s^2}, \end{equation}

where $z_j$ is one of the three solutions of $z_j^3 = 1$, and $s$ is one of the four solutions of

(3.27)\begin{equation} 25(x^2+y^2)s^4 + 20h^2 x s^3 + (4h^2+25x^2+20y^2)h^2 s^2 + 20h^4 x s + 4 h^4 (h^2+y^2) = 0. \end{equation}

This produces twenty-four potential late-order contributions, corresponding to the choice of sign and $z_j$ in (3.26) and the four solutions of (3.27). Twelve of these solutions are spurious, introduced by squaring both sides of an equation in the algebraic manipulations. This leaves twelve solutions, which appear as six complex conjugate pairs. Three of these pairs do not demonstrate Stokes switching, as they are exponentially large on the curve $\mathrm {Im}(\chi ) = 0$; they must therefore be inactive on the surface, and do not contribute to the free-surface behaviour. This leaves three singulant pairs that produce waves on the free surface. We will denote these as $\chi _j$ and $\overline {\chi }_j$ for $j = 1, 2, 3$, where the bar represents complex conjugation.

The behaviour of $\chi _1$ is illustrated in figure 9. The surface waves are absent directly downstream from the obstacle. The surface contains a Stokes curve that passes through the origin. This Stokes curve causes elastic waves to be switched on upstream from the obstacle. From direct algebraic computation, we find that $\chi \sim (1-\mathrm {i} x)/h$ in the limit that $x \rightarrow -\infty$. These elastic waves do not decay exponentially in space, although they will have algebraic spatial decay due to the prefactor, calculated below.

Figure 9. The (a) real and (b) imaginary parts of the singulant $\chi _1$, corresponding to algebraically decaying elastic waves. The Stokes curve, satisfying $\mathrm {Re}(\chi_1 ) > 0$ and $\mathrm {Im}(\chi _1) = 0$, is depicted in (b). The waves are not present in a region downstream from the obstacle, and appear as the Stokes curve is crossed into the region ahead of the obstacle. The equal phase lines in (b) illustrate the shape of the surface waves.

The behaviour of $\chi _2$ is illustrated in figure 10. The surface waves are also absent directly downstream from the obstacle. In fact, we see that this does not depend on the radiation condition; instead, there is an anti-Stokes curve on the surface. If the surface waves were present on this side of the Stokes curve, then they would become exponentially large on the downstream side of the anti-Stokes curve. We therefore see that the waves are present only on the downstream side of the Stokes curve. Importantly, $\mathrm {Re}(\chi )$ grows monotonically in the negative $x$-direction, becoming arbitrarily large as $x \rightarrow -\infty$. This means that the waves decay exponentially in space. We will therefore not consider these waves in the subsequent analysis.

Figure 10. The (a) real and (b) imaginary parts of the singulant $\chi _2$, corresponding to one set of exponentially decaying elastic waves. The anti-Stokes curve satisfies $\mathrm {Re}(\chi _2) = 0$, and is depicted as a dashed line. On the inside of this curve (where $\mathrm {Re}(\chi _2) < 0$), the exponential term would be large. Consequently, the remainder must be inactive in a region containing the anti-Stokes curve. The Stokes curve, satisfying $\mathrm {Re}(\chi _2) > 0$ and $\mathrm {Im}(\chi _2) = 0$, is depicted in (b). We see that the exponential must be switched on as this curve is crossed in a direction heading away from the origin. As the real part of $\chi _2$ increases without bound as $x \rightarrow -\infty$, these waves must decay exponentially in space, which is too rapid to have an observable physical effect. The behaviour of $\chi _3$ can be obtained by the mapping $x \mapsto -x$.

Finally, we find that $\chi _3$ is identical to $\chi _2$ reflected around the $y$-axis. This is consistent with the two-dimensional result, seen in Lustri et al. (2019), in which the surface contains symmetric waves emerging in both directions above the obstacle, which decay exponentially in space.

3.4.2. Calculating the prefactor

We follow an analysis similar to previous work on gravity and capillary waves in order to determine the prefactor expression $\varPhi$, and hence $\varXi$. These calculations are quite technical, and the details are included in § B.1. The prefactor $\varPhi$ is given by

(3.28)\begin{equation} \varPhi\!=\! \frac{s^{1/3}\sqrt{2}}{4{\rm \pi}^{3/2} h^{5/6}}\left[1\!-\!\frac{5 h^2 (2 h^2 \!-\! s^2) (s \!-\! x)}{3s (2 h^2 \!+\! 5 s^2) }\right]^{{({-}1)^{1/3} (4h^{10}-5h^8s^2-30h^2s^8-48s^{10})}/{30 h^{8/3}s^{16/3}(2h^2-s^2)}}, \end{equation}

where $s$ is the solution of (3.27) corresponding to the singulant illustrated in figure 9.

Finally, to find $\gamma$, we ensure that the strength of the singularity in the late-order behaviour $\phi ^{(n)}$ given in (3.17a,b) is consistent with the leading-order behaviour $\phi ^{(0)}$, which has strength $1/2$. It is clear from the recurrence relation (3.13) that the strength of the singularity will increase by three between $\phi ^{(n-1)}$ and $\phi ^{(n)}$. This implies that near the singularity at $x^2 + y^2 + h^2 = 0$,

(3.29)\begin{equation} \frac{\varPhi\,\varGamma(\gamma)}{\chi^{\gamma}} \rightarrow \frac{\alpha(x,y)}{(x^2+y^2+h^2)^{1/2}}, \end{equation}

where $\alpha$ is of order 1 in the limit. From (3.28), we see that the prefactor is also order 1 in this limit. A local analysis near the singularity (performed in (B13)) shows that $1/\chi$ contains a singularity with strength one at $x^2 + y^2 + h^2 = 0$. Matching the order of the expressions in (3.29) therefore gives $\gamma = 1/2$. We have therefore completely described the late-order terms in (3.17a,b), where (3.23a,b) is used to determine the value of $\varXi$, and hence the behaviour of the free-surface waves.

In § B.3, we use the late-order terms in (3.17a,b) to apply the matched asymptotic expansion methodology of Olde Daalhuis et al. (1995). We optimally truncate the asymptotic series and identify the Stokes curves. Finally, we use a matched asymptotic expansion analysis on the truncation remainder to compute the exponentially small contribution to the free-surface behaviour that appears across the Stokes lines. This analysis in § B.3 follows steps similar to the equivalent analysis in Lustri & Chapman (2013) and Lustri et al. (2019).

Using this method, we find that the exponentially small contributions to the fluid potential (denoted $\phi _{{exp}}$) and free-surface position (denoted $\xi _{{exp}}$) as $\epsilon \to 0$ are switched in the region to the left of the Stokes curve shown in figure 9. In the region where the exponentially small contributions are present, they are given by

(3.30a,b)\begin{equation} \phi_{{exp}} \sim \frac{2{\rm \pi}\mathrm{i} \varPhi}{3\sqrt{\epsilon}}\,\mathrm{e}^{-\chi_1/\epsilon} + \mathrm{c.c.},\quad \xi_{{exp}} \sim \frac{2{\rm \pi}\mathrm{i} \varXi}{3\sqrt{\epsilon}}\,\mathrm{e}^{-\chi_1/\epsilon} + \mathrm{c.c.},\end{equation}

where c.c. denotes the complex conjugate contribution. In particular, the expression for $\xi _{{exp}}$ contains exponentially small oscillations as $\epsilon \to 0$, representing the elastic ripples on the free surface.

We note that (3.30b) is a leading-order expression for the exponentially small waves, with algebraic corrections to the prefactors $\varPhi$ and $\varXi$ omitted. The solution also contains contributions from $\chi _2$ and $\chi _3$ with the same form, but these decay exponentially in space, and are therefore exponentially smaller in amplitude as $\epsilon \to 0$ than the elastic waves caused by $\chi _1$, as well as the neglected correction terms. We therefore omit these contributions from the asymptotic expression.

3.5. Results and comparison

Along the curve $y=0$ for $x < 0$, we have $s = \mathrm {i} h$ and $\chi =h + \mathrm {i} x$. We therefore evaluate the free-surface position to be

(3.31)\begin{equation} \xi_{{exp}} \sim{-}\frac{ ({-}1)^{1/6} 3^{3/10}}{3 h^{1/5} \sqrt{2{\rm \pi}\epsilon}\,(8h + 5 \mathrm{i} x)^{3/10}}\,\mathrm{e}^{-(h + \mathrm{i} x)/\epsilon} + \textrm{c.c.} \quad \mathrm{as} \ \epsilon \to -\infty, \end{equation}

where c.c. denotes the complex conjugate contribution. In the limit that $x$ becomes large and negative, we find that the amplitude of the waves on $y = 0$ is is given by

(3.32)\begin{equation} \mathrm{amplitude} \sim \frac{1}{3 h^{1/5}}\sqrt{\frac{2}{{\rm \pi}\epsilon}} \left(\frac{3}{5|x|}\right)^{3/10}\mathrm{e}^{{-}h/\epsilon} \quad \mathrm{as} \ x \to -\infty,\, \epsilon \rightarrow 0.\end{equation}

This provides us with a quantity that we may use to check the accuracy of the asymptotic approximation. We compare the amplitude of the asymptotic results with those of numerically calculated free-surface profiles obtained by solving the linearized system (3.6)–(3.9). These results were obtained using an adaptation of the method described in Lustri & Chapman (2013) and Lustri et al. (2019), which consists of expressing the free-surface behaviour in terms of Fourier inversion integrals, and evaluating the double integral numerically on a fixed domain.

In figure 11, we illustrate the scaled numerical amplitude (circles) against the asymptotic prediction from (3.32), computed for $h = 1$ over a range of $\epsilon$ values. The amplitude is scaled by $|x|^{3/10}$, so that it tends to a constant as $x \to -\infty$. The numerical amplitude is taken by determining the scaled amplitude for sufficiently large negative values of $x$ that the scaled amplitude does not display significant variation. There is strong agreement between the asymptotic predictions and numerical results. For values of $\epsilon$ smaller than those depicted, it become numerically challenging to compute the wave behaviour, due to the very small amplitude of the resulting waves.

Figure 11. Comparison between asymptotic predictions of the amplitude using (3.32) and numerical calculations, for the case $h=1$. The amplitude is scaled by $|x|^{3/10}$ so that it tends to a constant value as $x \rightarrow -\infty$.

4. Discussion and conclusions

In this paper, we studied behaviour of the waves that form on an elastic sheet resting on an inviscid flow stream containing a submerged obstacle in several different regimes. In each regime that we considered, the surface waves are exponentially small in the limit of small bending stiffness, and we therefore required exponential asymptotic techniques in order to study the wave behaviour. As these waves share many similarities with surface-tension-driven capillary waves, we used techniques applied in Trinh & Chapman (2013a,b) for two-dimensional gravity-capillary waves, and Lustri et al. (2019) for three-dimensional capillary waves. Using exponential asymptotics, we identified the regions in which the surface waves formed, and calculated a mathematical expression for the wave behaviour. In each flow configuration, we determined that the surface waves appeared as Stokes curves were crossed, and by understanding the behaviour of these Stokes curves, it was possible to classify the types of waves that could appear on the elastic sheet.

We first studied the behaviour of flexural-gravity waves on linearized flow over a small step in two dimensions. In this regime, the elastic sheet behaviour depended on a particular parameter $\beta _c(\tau )$, which related the relative sizes of the Froude number and bending stiffness parameter. If the parameter $\beta$ is less than this critical value, then the four different wave contributions on the elastic sheet decay away from the obstacle, meaning that both the upstream and downstream regions do not contain any waves. If $\beta$ exceeds this critical value, then constant-amplitude elastic waves propagate upstream from the obstacle, and constant-amplitude gravity waves propagate downstream. In this regime, the sheet behaviour also contains two decaying wave contributions, but these are small compared to the constant-amplitude waves.

We did not perform a full exponential asymptotic analysis of the two-dimensional nonlinear problem, but instead used the singulant behaviour and Stokes structure of the problem to predict the types of wave that could form on the elastic sheet. Flow over a submerged step in this case contains a third intermediate regime, in addition to the two regimes from the linearized problem. For flow over an upward step, there are two distinct critical values: one critical value of $\beta$ that determines whether constant-amplitude downstream gravity waves are present, and a larger critical value of $\beta$ that determines whether constant-amplitude upstream elastic waves are present. It is therefore possible to construct flow geometries that contain constant-amplitude downstream gravity waves, but where all upstream waves decay in space. The converse is true for a downward step: it is possible to construct geometries that contain constant-amplitude upstream elastic waves, but all downstream waves decay in space.

In the two-dimensional geometries considered in this study, we considered the behaviour of systems in which gravitational and elastic effects have similar strength, while neglecting surface tension effects. It would also be interesting to calculate the wave behaviour while incorporating interactions between surface tension of the fluid and bending effects. This interaction was studied experimentally in Deike, Bacri & Falcon (2013) and Deike, Berhanu & Falcon (2017), and a similar form has been used to study compressive effects on elastic sheets in Das, Sahoo & Meylan (2018). In this case, both capillary and elastic waves would be expected to propagate upstream, and the behaviour in the elastic sheet would be more straightforward to validate computationally.

Finally, we studied hydroelastic waves that form on linearized flow in three dimensions over a submerged obstacle. We considered a regime in which the wave behaviour is caused solely by the elastic restoring force. Using exponential asymptotics, we calculated the wave behaviour, and found that it is consistent with the two-dimensional case from Lustri et al. (2020). The flow contained four wave contributions, two of which decay exponentially in space, and one of which produced visible upstream elastic waves that decay algebraically. We validated these calculations against numerical results. The wave patterns are visible similar to the capillary waves from Lustri et al. (2019), although the algebraic decay rate of the wave amplitude differs between elastic and capillary waves, and show qualitative agreement with the experimental results of Ono-dit-Biot et al. (2019).

It is natural to consider whether this analysis can be extended to flexural-gravity waves in three dimensions. In this case, we would replace (3.3) with

(4.1)\begin{equation} \frac{\beta\epsilon}{2}\left(|\boldsymbol{\nabla}\phi|^2-1\right) +\xi + \beta\tau\epsilon^4(\xi_{xxxx} + 2 \xi_{xxyy} + \xi_{yyyy}) = 0,\quad z = \xi(x,y), \end{equation}

where $F^2 = \beta \epsilon$ and $\sigma = \beta \tau \epsilon ^4$, as in the two-dimensional problem. For the linearized problem, repeating the late-order analysis in a similar fashion gives

(4.2)\begin{equation} \beta \chi_x^4 + (\chi_x^2 + \chi_y)^2[1 + \beta\tau(\chi_x^2 + \chi_y^2)^2]^2 = 0, \end{equation}

with the boundary condition $\chi = 0$ on $x^2 + y^2 + h^2 = 0$. This equation is challenging to solve using analytical methods, as there are eight valid solution sheets in the analytically continued free surface, which has $x\in \mathbb {C}$ and $y\in \mathbb {C}$. Solving this equation directly would likely require more complicated ray-tracing methods such as those developed for nonlinear three-dimensional flow in Johnson-Llambias (2022), and would therefore be unlikely to identify convenient closed-form solutions such as (3.30a,b). Nonetheless, these computations would likely still be useful, as exponential asymptotic methods are a convenient method for isolating particular wave behaviours and studying the waves directly.