2.1. Governing equations

The fluid flow is governed by the (two-dimensional) incompressible Navier–Stokes equations, where ${\boldsymbol u}=(u,v)$ and $p$ denote the dimensionless velocity and pressure of the fluid, respectively, in the form

(2.2a)\begin{gather} {u}_x+{v}_y=0, \end{gather}

(2.2b)\begin{gather}{u}_t+u{u}_x+v{u}_y={-}{p}_x+{R}^{{-}1}(u_{xx} + u_{yy}), \end{gather}

(2.2c)\begin{gather}{v}_t+u{v}_x+v{v}_y={-}{p}_y+{R}^{{-}1}(v_{xx} + v_{yy}). \end{gather}

On the rigid wall we apply boundary conditions of no slip and no penetration in the form

(2.2d)\begin{equation} u=0, \quad v=0, \quad (y=0),\end{equation}

while on the compliant wall we apply kinematic and no-slip conditions in the form

(2.2e)\begin{equation} u={-}\frac{h_t h_x}{1+h_x^2}, \quad v=\frac{h_t}{1+h_x^2}, \quad (y=h).\end{equation}

The shape of the flexible wall is governed by a normal stress balance, in the form

(2.2f)\begin{align} &p-\frac{2}{R(1+h_x^2)}(v_y-h_x(u_y+v_x)+h_x^2u_x) \nonumber\\ &\quad=p_e(x)-\frac{Th_{xx}}{(1+h_x^2)^{{3}/{2}}}+ \frac{mh_{tt}}{({1+h_x^2})^{1/2}}+\frac{dh_t}{({1+h_x^2})^{1/2}}, \quad (y=h). \end{align}

This wall model includes the same features as Benjamin (Reference Benjamin1960, Reference Benjamin1963) and Landahl (Reference Landahl1962). We note that this balance includes the normal viscous stress (terms appearing with coefficient $R^{-1}$).

2.2. Nonlinear energy equation

To analyse the energy budget of the modes of instability in this collapsible channel system we compute the energy budget by taking the scalar product of the momentum equations (2.2b,c) with the fluid velocity $\boldsymbol {u}$ in the form

(2.3a)\begin{align} &(u{u}_t+v{v}_t)+u({u}{u}_x+v{v}_x)+v (u{u}_y+{v}{v}_y)+\tfrac{1}{2}(u^2+v^2) (u_x+v_y) \nonumber\\ &\quad ={-}u{p}_x-v{p}_y+R^{{-}1}(uu_{xx} + uu_{yy}+vv_{xx} + vv_{yy} )-p(u_x+v_y). \end{align}

This energy equation can be manipulated into the familiar form

(2.3b)\begin{equation} {\mathcal{F}}+{\mathcal{P}}={\mathcal{K}}+{\mathcal{E}}+{\mathcal{D}},\end{equation}

where

(2.3c)\begin{gather} {\mathcal{F}}={-}\left(\frac{1}{2}\int_{0}^{h}u (u^2+v^2)\,{\rm d} y\right)_x, \end{gather}

(2.3d)\begin{gather}{\mathcal{P}}={-}\left(\int_{0}^{h} up\,{\rm d} y\right)_x, \end{gather}

(2.3e)\begin{gather}{\mathcal{K}}=\left(\frac{1}{2}\int_{0}^{h}(u^2+v^2)\,{\rm d} y\right)_t, \end{gather}

(2.3f)\begin{align} {\mathcal{E}}&={-}{\mathcal{W}}=h_t \left[p-\frac{2}{R(1+{h_x}^2)} (v_y-h_x(u_y+v_x)+{h_x^2} u_x)\right]^{y=h}={-} [h_t({\boldsymbol{\sigma}}{\boldsymbol n})\boldsymbol{\cdot} {\boldsymbol n}]^{y=h} \nonumber\\ &=p_e h_t - T\frac{h_t h_{xx}}{(1+h_x^2)^{{3}/{2}}}+m\frac{h_t h_{tt}}{({1+h_x^2})^{1/2}}+d\frac{h_t^2}{({1 +h_x^2})^{1/2}}, \end{align}

(2.3g)\begin{equation} {\mathcal{D}}={R}^{{-}1}\left\{\int_0^h(2u_x^2+2v_y^2+(u_y+v_x)^2)\,{\rm d} y-\left(\int_0^h(2uu_x+vv_x+vu_y)\,{\rm d} y\right)_x\right\}, \end{equation}

where ${\boldsymbol {\sigma }}$ is the Newtonian stress tensor and ${\boldsymbol n}$ is the outward pointing unit normal to the flexible wall (see figure 1). Here, ${\mathcal {F}}$ is the kinetic energy flux extracted from the mean flow, $\mathcal {P}$ is the rate of working of the axial pressure gradient, ${\mathcal {K}}$ is the rate of change of kinetic energy, ${\mathcal {W}}=-{\mathcal {E}}$ is the rate of working of fluid normal stress on the wall and ${\mathcal {D}}$ is the rate of energy loss due to fluid viscosity in the bulk. Note that we re-write ${\mathcal {E}}$ by substituting the normal stress balance (2.2f) and the boundary conditions (2.2e), which then expresses the work done by the fluid on the wall in terms of the wall parameters. Note also that previous studies (e.g. Carpenter & Gajjar Reference Carpenter and Gajjar1990; Huang Reference Huang1998) have considered the work done by fluid pressure on the wall, whereas in this case we consider the work done by the normal stress, including the viscous contributions (see also Lebbal et al. Reference Lebbal, Alizard and Pier2022).

2.3. Linearisation around the uniform basic state

We assume the velocity and pressure of the steady flow when the elastic wall is flat to be of the familiar (fully developed) Poiseuille form

(2.4)\begin{equation} (u,v)=(U(y),0)=(6y(1-y),0),\quad p=p_e=P(x)=C-\frac{12}{R}x, \quad h=1, \end{equation}

where $C$ is a constant, which is a solution to (2.2) for all $C$. We note that the external pressure is not spatially uniform; we are applying an external pressure gradient to maintain a flat base state. Such an assumption was previously applied by Stewart et al. (Reference Stewart, Waters and Jensen2010b) and is consistent with earlier work which ignored the pressure gradient associated with the mean flow in the normal stress balance (e.g. Carpenter & Garrad Reference Carpenter and Garrad1985; Davies & Carpenter Reference Davies and Carpenter1997a).

We add a small perturbation to this steady flow in the form

(2.5)\begin{equation} (u,v,p,h)=(U(y),0,P(x),1)+\varepsilon (\hat{\phi}_y,-\hat\phi _x,\hat{p},\hat{\eta})+ \varepsilon^2(\breve\phi _y,-\breve\phi _x,\breve{p},\breve{\eta}) +O(\varepsilon^3),\end{equation}

where $\varepsilon \ll 1$ is the perturbation amplitude, and $\hat \phi$ and $\breve \phi$ are the perturbation streamfunctions at $O(\varepsilon )$ and $O(\varepsilon ^2)$, respectively. Note that conservation of mass (2.2a) is satisfied automatically at each order.

We substitute the (2.5) into the Navier–Stokes equations (2.2b,c), the boundary conditions (2.2d,e) and the normal stress balance (2.2f) and examine successive orders in the perturbation amplitude $\varepsilon$.

At $O(1)$ these equations are satisfied automatically, while at $O(\varepsilon )$ the governing equations of the system become

(2.6a)\begin{gather} \hat\phi _{yt}+U\hat\phi _{xy}-\hat\phi _{x}U_y+\hat p_x-R^{{-}1}(\hat\phi _{xxy}+\hat\phi _{yyy})=0, \end{gather}

(2.6b)\begin{gather}\hat\phi _{xt}+U\hat\phi _{xx}-\hat p_y-R^{{-}1}(\hat\phi _{xxx}+\hat\phi _{xyy})=0, \end{gather}

subject to

(2.6c)\begin{gather} \hat\phi_y=0, \quad \hat\phi_x=0, \quad (y=0), \end{gather}

(2.6d)\begin{gather}\hat\phi_y+U_y(1)\hat\eta=0,\quad \hat\phi_x + \hat\eta_t=0,\quad (y=1), \end{gather}

while the normal stress balance takes the form

(2.6e)\begin{equation} \hat p +{2}{R}^{{-}1}(\hat\phi_{y} + U_y(1)\hat\eta)_x + T\hat\eta_{xx}-m\hat\eta_{tt}-d\hat \eta_t =0, \quad (y=1).\end{equation}

Here, we have explicitly included the normal viscous stress contributions in (2.6e), although these cancel identically at this order due to the no-slip boundary condition (2.6d).

Similarly, at $O(\varepsilon ^2)$, the governing equations of the system become

(2.7a)\begin{gather} \breve\phi_{yt}+U\breve\phi_{xy}-U_y\breve\phi_x+\breve p_x-R^{{-}1}(\breve\phi_{xxy}+\breve\phi_{yyy}) =\hat\phi_x\hat\phi_{yy}-\hat\phi_y\hat\phi_{xy}, \end{gather}

(2.7b)\begin{gather}\breve\phi_{xt}+U\breve\phi_{xx}-\breve{p}_y-R^{{-}1} (\breve\phi_{xxx}+\breve\phi_{xyy}) =\hat\phi_x\hat{\phi}_{xy}-\hat\phi_y\hat\phi_{xx}, \end{gather}

subject to boundary conditions

(2.7c)\begin{gather} \breve\phi_y=0, \quad \breve\phi_x=0, \quad (y=0), \end{gather}

(2.7d)\begin{gather}\breve\phi_y + U_y(1)\breve\eta ={-}\hat\phi_{yy}\hat\eta - \tfrac{1}{2}U_{yy}(1) {\hat \eta}^2 - \hat\eta_x\hat\eta_t,\quad \breve\phi_x+\breve\eta_t={-}\hat \eta\hat\phi_{xy}, \quad (y=1), \end{gather}

and the normal stress balance

(2.7e)\begin{align} &\breve p+{2}{R^{{-}1}}(\breve\phi_{xy}+U_y(1)\breve\eta_x) +T\breve\eta_{xx}-m\breve\eta_{tt}-{\rm d} \breve \eta_t \nonumber\\ &\quad ={-}\hat\eta\hat p_y-{2}{R^{{-}1}}(\hat\phi_{xyy}\hat\eta +U_{yy}(1)\hat\eta\hat\eta_x+(\hat\phi_{yy}-\hat\phi_{xx}) \hat\eta_x),\quad (y=1). \end{align}

2.4. The first-order equations

To solve the first-order governing equations, we assume that all hatted variables are wave like in the form

(2.8)\begin{equation} {\hat f}=\hat{\hat{f}}(y)\, {\rm e}^{{\rm i}(kx-\omega t)}+(\,\hat{\hat{f}}(y)\, {\rm e}^{{\rm i}(kx-\omega t)})^*,\end{equation}

where $k$ is the wavenumber and $\omega$ is the corresponding frequency of the perturbation, while $^*$ denotes the complex conjugate. We also express our neutral stability curves in terms of the perturbation wavespeed, $c=\omega /k$ for $k$ real (Whitham Reference Whitham1974).

Substituting (2.8) into the first-order governing equations (2.6a)–(2.6e), we have

(2.9a)\begin{gather} \left(U-\frac{\omega}{k}\right)\hat{\hat{\phi}}_y- \hat{\hat{\phi}}U_y={-}\hat{\hat{p}}+\frac{1}{ikR} ({-}k^2\hat{\hat{\phi}}_y+\hat{\hat{\phi}}_{yyy}), \end{gather}

(2.9b)\begin{gather}k^2\left(U-\frac{\omega}{k}\right)\hat{\hat{\phi}} ={-}\hat{\hat{p}}_y+\frac{1}{ikR}({-}k^4\hat{\hat{\phi}} +k^2\hat{\hat{\phi}}_{yy}), \end{gather}

subject to boundary conditions

(2.9c)\begin{gather} \hat{\hat{\phi}}=0, \quad \hat{\hat{\phi}}_y=0,\quad (y=0), \end{gather}

(2.9d)\begin{gather}\hat{\hat{\phi}}=\frac{\omega}{k}\hat{\hat{\eta}}, \quad \hat{\hat{\phi}}_y={-}U_y(1)\hat{\hat{\eta}},\quad (y=1), \end{gather}

and the normal stress balance

(2.9e)\begin{equation} \hat{\hat{p}}=\frac{1}{ikR}({-}k^2\hat{\hat{\phi}}_y+ \hat{\hat{\phi}}_{yyy})=(Tk^2-m\omega^2-id\omega) \hat{\hat{\eta}}, \quad (y=1).\end{equation}

Combining (2.9a) and (2.9b) by eliminating $\hat {\hat {p}}$, we obtain the classical Orr–Sommerfeld equation (see below, Drazin & Reid Reference Drazin and Reid1981; Schmid & Henningson Reference Schmid and Henningson2001).

We adopt two complementary approaches to solve this linear system. First, we use the no-slip boundary condition (2.9d) to eliminate the wall amplitude $\hat {\hat {\eta }}$ to derive a (closed) fourth-order homogeneous system for ${\hat {\hat \phi }}(y)$ in the form

(2.10a)\begin{equation} \left(U-\frac{\omega}{k}\right) (\hat{\hat{\phi}}_{yy}-k^2\hat{\hat{\phi}})- U_{yy}\hat{\hat{\phi}}-\frac{1}{ikR} (\hat{\hat{\phi}}_{yyyy}-2k^2\hat{\hat{\phi}}_{yy} +k^4\hat{\hat{\phi}})=0, \end{equation}

subject to four boundary conditions

(2.10b)\begin{gather} \hat{\hat{\phi}}(0)=0, \quad \hat{\hat{\phi}}_y(0)=0, \end{gather}

(2.10c)\begin{gather}\omega\hat{\hat{\phi}}_y(1) + U_y(1)k \hat{\hat{\phi}}(1)=0, \end{gather}

(2.10d)\begin{gather}\frac{1}{i k R}U_y(1)({-}k^2\hat{\hat{\phi}}_y(1)+\hat{\hat{\phi}}_{yyy}(1) ) + (Tk^2 - m\omega^2-id\omega)\hat{\hat{\phi}}_y(1)=0. \end{gather}

For fixed Reynolds number, $R$, and wavenumber, $k$, we solve the eigenvalue problem numerically by discretising (2.10) using a Chebyshev spectral method similar to that used by Stewart et al. (Reference Stewart, Waters and Jensen2010b) modified to include the wall mass. For a given value of the wavenumber and the other model parameters ($R$, $T$, $m$ and $d$), for $N$ Chebyshev modes the numerical method isolates eigenvalues $\omega _i$ and corresponding eigenfunctions ${\hat {\hat \phi }}_i(y;k,\omega )$ for $i=1,2,\ldots,N$. Typical eigenvalue spectra are shown in figure 6(c,d), where the eigenvalues can be divided into the familiar A, P and S branches (Mack Reference Mack1976; Drazin & Reid Reference Drazin and Reid1981; Schmid & Henningson Reference Schmid and Henningson2001). As the number of Chebyshev modes increases the additional eigenvalues appear with increasingly large and negative imaginary part (along the S branch) and the portion of the spectrum close to neutral stability becomes well resolved. See Wang (Reference Wang2019) for further details. Knowledge of the entire eigenvalue spectrum (as opposed to just neutrally stable modes) serves to elucidate the mode interactions discussed in later sections. For the plots in this paper we normalise the eigenfunction so that $\hat {\hat {\phi }}_y(1)=-U_y(1)=6$, consistent with the inhomogeneous formulation (2.11). We focus particular attention on neutrally stable solutions, where ${\rm Im} (\omega )=0$, which we isolate using a bisection method and trace various slices through the parameter space. In particular, we compute the neutrally stable wavenumber, denoted $k_n$, corresponding (real) neutrally stable frequency, denoted $\omega _n$ and corresponding (real) wavespeed, denoted $c_n=\omega _n/k_n$, as a function of the model parameters in figures 2–10.

To make analytical progress we also adopt an alternative approach which separates computation of the eigenfunction from the stability calculation. In this case we first normalise the eigenfunction based on the amplitude of the wall profile, scaling the perturbation streamfunction according to $\hat {\hat {\phi }} = \hat {\hat {\eta }}{\tilde {\hat \phi }}$ and express (2.9) as a closed inhomogeneous system for ${\tilde {\hat \phi }}(y)$ alone in the form

(2.11a)\begin{equation} \left(U-\frac{\omega}{k}\right)(\tilde{\hat{\phi}}_{yy} -k^2\tilde{\hat{\phi}})-U_{yy}\tilde{\hat{\phi}}- \frac{1}{ikR}(\tilde{\hat{\phi}}_{yyyy}-2k^2 \tilde{\hat{\phi}}_{yy}+k^4\tilde{\hat{\phi}})=0, \end{equation}

subject to

(2.11b)\begin{gather} \tilde{\hat{\phi}}(0)=0, \quad \tilde{\hat{\phi}}_y(0)=0, \end{gather}

(2.11c)\begin{gather}\tilde{\hat{\phi}}(1)=\frac{\omega}{k}, \quad \tilde{\hat{\phi}}_y(1)={-}U_y(1). \end{gather}

In this way, the streamfunction ${\tilde {\hat \phi }}(y)$ can be computed for any $k$ and $\omega$, independent of the stability criteria. Landahl (Reference Landahl1962) achieved an analogous separation in his (more approximate) approach by expressing the problem in terms of the wall admittance, which facilitated calculation of the activation energy. Given a solution of (2.11), following Cairns (Reference Cairns1979), the normal stress balance on the flexible wall can be expressed as a dispersion relation of the system in the form

(2.11d)\begin{equation} D(\omega,k)\equiv\frac{1}{ikR}(\tilde{\hat{\phi}}_{yyy}(1)- k^2\tilde{\hat{\phi}}_y(1))-(Tk^2-m\omega^2-id\omega),\end{equation}

which is zero for normal modes. Crucially, this approach formulates the dispersion relation as a scalar-valued function involving the frequency $\omega$ and the wavenumber $k$ (in this case in terms of the normal stress balance on the wall), rather than as a inhomogeneous eigenvalue problem which must be solved numerically (e.g. Davies & Carpenter Reference Davies and Carpenter1997a; Stewart et al. Reference Stewart, Waters and Jensen2010b); this scalar dispersion relation has two components: the first group of terms arises due to the viscous normal stress on the wall, while the second group of terms arises from the contributions to the Euler–Bernoulli beam. The inviscid dispersion relation (see Appendix A) does not have the viscous terms but contains an extra term due to inertial contribution to the pressure, which cancels identically here due to the no-slip condition on the flexible wall. We use this dispersion relation to compute the activation energy of the neutrally stable surface-based modes following Landahl (Reference Landahl1962) and Cairns (Reference Cairns1979).

Note that our choice of normalisation for the numerical solutions of the homogeneous system (2.10) means that at neutral stability the two streamfunction profiles are identical ($\hat {\hat {\phi }}(y;\omega _n,k_n)=\tilde {\hat {\phi }}(y;\omega _n,k_n)$). For notational simplicity, we henceforth denote $\tilde {\hat {\phi }}(y)=\phi (y)$.

2.5. The second-order equations

At second order we express the variables in the form

(2.12)\begin{equation} \breve g(x,y,t)=\breve{\breve{g}}(y)\, {\rm e}^{2{\rm i}(kx-\omega t)}+\bar{{\breve{g}}}(y)+(\breve{\breve{g}}(y)\, {\rm e}^{2{\rm i}(kx-\omega t)})^*. \end{equation}

Substituting (2.12) into the second-order equations (2.7), we consider steady streaming terms which are independent of both $x$ and $t$ (denoted with an additional overline) to obtain a system which depends only on $y$, in the form

(2.13a)\begin{equation} {-\bar{\breve p}_x +} \frac{1}{R}\bar{\breve{\phi}}_{yyy} = \overline{\hat{\phi}_y \hat{\phi}_{xy}} {- \overline{\hat{\phi}_x \hat{\phi}_{yy}}}, \quad \bar{{\breve p}}_y=0, \quad (0\le y \le 1),\end{equation}

subject to

(2.13b)\begin{equation} \bar{\breve{\phi}}_{y}(0)=0, \quad \bar{\breve{\phi}}_{y}(1)={-} U_{y}(1)\bar{{\breve \eta}} - \tfrac{1}{2}U_{yy}(1)\overline{{\hat \eta}^2} - \overline{{\hat \phi}_{yy}(1) {\hat \eta}} - \overline{{\hat \eta}_t {\hat \eta}_x},\end{equation}

noting that $x$ and $t$ independent terms occur at this order from products of the first-order eigenfunctions. Note also that the term $\bar {\breve p}_x$ is independent of $x$, representing the steady modification to the pressure gradient along the channel. The corresponding kinematic boundary condition is identically zero. Steady streaming can either result in a modification to the steady pressure gradient ($\bar {\breve p}_x$), a modification to the steady flux along the channel ($\bar {\breve q}=\bar {\breve \phi }(1)-\bar {\breve \phi }(0)$) or the mean wall position (${\bar {\breve \eta }}$), only one of which can be fixed by solving the equations. In this case we are assuming a fixed baseline flow rate which imposes $\bar {\breve q}=0$ and ${\bar {\breve \eta }}=0$ to ensure the total mass of fluid in the channel remains fixed. By integrating (2.13a) and applying boundary conditions (2.13b) we obtain a prediction of the streamfunction profile $\bar {\breve {\phi }}_{y}$ and pressure gradient $\overline {\breve {p}}_x$ in terms the numerically computed eigenfunctions from the first-order problem (see Wang (Reference Wang2019), for details). This profile is used to calculate the full expression for the work done on the wall by the fluid in the time- and space-averaged energy budget (see § 2.6).

2.6. The Reynolds–Orr energy equation

We obtain the linearised energy equation by substituting the expansion (2.5) into the energy equation (2.3). The first non-trivial balance emerges at second order in perturbation amplitude, where we have

(2.14)\begin{align} &U\breve\phi_{yt}+U^2\breve\phi_{xy}-UU_y\breve{\phi}_x +U\breve{p}_x+P_x\breve{\phi}_y-{R}^{{-}1}(U\breve{\phi}_{yxx} +U\breve{\phi}_{yyy}+U_{yy}\breve{\phi}_y)\nonumber\\ &\quad+(\hat{\phi}_x\hat\phi_{xt}+\hat{\phi}_y\hat\phi_{yt}) + (U\hat{\phi}_x\hat\phi_{xx}+2U\hat{\phi}_y\hat\phi_{xy}-U_y \hat{\phi}_x\hat{\phi}_y+U\hat{\phi}_x\hat\phi_{yy}) \nonumber\\ &\quad+\hat{p}_x\hat{\phi}_y-\hat{p}_y\hat{\phi}_x-{R}^{{-}1} (\hat{\phi}_y\hat{\phi}_{yxx}+\hat\phi_y\hat\phi_{yyy}+ \hat\phi_x\hat\phi_{xxx}+\hat\phi_x\hat{\phi}_{xyy})=0. \end{align}

Domaradzki & Metcalfe (Reference Domaradzki and Metcalfe1987) (see also Stewart et al. Reference Stewart, Waters and Jensen2010b; Lebbal et al. Reference Lebbal, Alizard and Pier2022) demonstrated that the terms involving breved variables can be eliminated in favour of terms involving the first-order variables alone (using the $O(\epsilon ^2)$ governing equations), analogous to the derivation of the classical Reynolds–Orr equation (see details in Schmid & Henningson Reference Schmid and Henningson2001). The corresponding energy budget is obtained by integrating across the static channel $(0\le y\le 1)$ to obtain the work done by nonlinear Reynolds stresses, which was termed ${\hat {\mathcal {S}}}$ by Stewart et al. (Reference Stewart, Waters and Jensen2010b). However, in this approach, it is not then possible to explicitly construct the full expression for the work done by the fluid on the wall as some contributions are embedded within ${\hat {\mathcal {S}}}$. In this study, we adopt an alternative approach where we retain the full expressions for the terms in the energy equation (2.14) involving both hatted and breved variables. In this case the energy budget at $O(\varepsilon ^2)$ can be expressed as

(2.15a)\begin{equation} {\hat{\mathcal{F}}} + {\hat{\mathcal{P}}} ={\hat{\mathcal{K}}} + {\hat{\mathcal{E}}} + {\hat{\mathcal{D}}}, \end{equation}

where

(2.15b)\begin{gather} {\hat{\mathcal{F}}}=\left( -\int_{0}^{1}\frac{1}{2}U(3{\hat\phi_y}^2+{\hat\phi_x}^2)\,{\rm d} y-\int_0^1\frac{3}{2}\breve\phi_yU^2\,{\rm d} y\right)_x, \end{gather}

(2.15c)\begin{gather}{\hat{\mathcal{P}}}= \left(-\int_0^1\hat\phi_y\hat p\,{\rm d} y\right)_x- \int_0^1(P_x\breve\phi_y+{U{\breve p}_x})\,{\rm d} y, \end{gather}

(2.15d)\begin{gather}{\hat{\mathcal{K}}}= \left(\int_{0}^{1}\frac{1}{2} ({\hat\phi_x}^2+{\hat\phi_y}^2)\,{\rm d} y + \int_0^1 U\breve\phi_y\,{\rm d} y\right)_t, \end{gather}

(2.15e)\begin{align} {\hat{\mathcal{E}}}&={-}{\hat{\mathcal{W}}}={-}(T\hat\eta_t\hat\eta_x)_x+\tfrac{1}{2} (T\hat\eta_x^2+m\hat\eta_t^2)_t+d\hat\eta_t^2\nonumber\\ &\quad {-}{R}^{{-}1}[-{U_{yy}}(\hat\eta\hat\phi_x)_x + \hat\phi_y\hat\phi_{yy}-3\hat\phi_y\hat\phi_{xx} + U_y\breve\phi_y]^{y=1}, \end{align}

(2.15f)\begin{align} {\hat{\mathcal{D}}}&={R}^{{-}1}\int_0^1(4\hat\phi_{xy}^2+ (\hat\phi_{yy}-\hat\phi_{xx})^2+2U_y\breve\phi_{yy})\,{\rm d} y \nonumber\\ &\quad -{R}^{{-}1}\left(\int_0^1(2\hat\phi_y\hat\phi_{xy}+ \hat\phi_x\hat\phi_{xx}-\hat\phi_x\hat\phi_{yy} + 2U\breve\phi_{xy}+U_y\breve\phi_x)\,{\rm d} y\right)_x; \end{align}

where ${\hat {\mathcal {F}}}$ represents the net energy extracted from the mean flow across the channel, ${\hat {\mathcal {P}}}$ represents the work done by axial pressure forces over a wavelength, ${\hat {\mathcal {W}}}=-{\hat {\mathcal {E}}}$ represents the total work done by the fluid normal stress on the wall (note the wall elasticity parameters appear as the normal stress condition (2.6e) has been used to eliminate the fluid pressure on the wall) and ${\hat {\mathcal {D}}}$ represents the total viscous dissipation across the bulk of the channel.

Following Landahl (Reference Landahl1962), we average (2.15) over a wavelength of perturbation where it reduces to

(2.16)\begin{align} &\left(\int_{0}^{1}\frac{1}{2}({\hat\phi_x}^2+{\hat\phi_y}^2) \,{\rm d} y + \frac{1}{2}(T\hat\eta_x^2+m\hat\eta_t^2)\right)_t ={-}d(\hat\eta_t)^2 \nonumber\\ &\quad -{R}^{{-}1}\int_0^1(4\hat\phi_{xy}^2+ (\hat\phi_{yy}-\hat\phi_{xx})^2+2U_y\breve\phi_{yy}) \,{\rm d} y \nonumber\\ &\quad -{R}^{{-}1}[\hat\phi_y\hat\phi_{yy}-3\hat\phi_y\hat\phi_{xx} + U_y\breve\phi_y]^{y=1} ; \end{align}

the term on the left-hand side is the rate of change of the activation energy of the system (sum of kinetic and elastic energies for the fluid and the wall), while the terms on the right-hand side represent (non-conservative) energy losses due to the wall damping, fluid viscosity across the bulk and viscous stresses exerted on the flexible wall, respectively. In a conservative system (with no wall damping or viscous dissipation) the activation energy is constant and must always be positive definite in this case. This suggests that undamped modes must always be class B (according to Benjamin's framework) in the limit of large Reynolds number. However, in the fully inviscid limit (Appendix A) the activation energy can be negative due to the extra (inertial) contribution to the fluid pressure on the wall which cancels in the viscous problem due to the no-slip condition. However, this approach cannot be used as a measure of activation energy for neutrally stable modes in a viscous system, since this formulation neglects the contribution of viscosity in establishing the neutrally stable mode. We propose an alternative method to construct the sign of the activation energy of neutrally stable viscous modes in the following subsection.

Instead, taking the average of (2.15) over both a period and wavelength (again denoted with an overbar), the averaged energy budget (2.15) reduces to

(2.17a)\begin{equation} \bar{\hat{\mathcal{P}}} =\bar{\hat{\mathcal{E}}} + \bar{\hat{\mathcal{D}}}, \end{equation}

where

(2.17b)\begin{gather} {\bar{\hat{\mathcal{P}}}}={-}\int_0^1 U \overline{\breve{p}}_x\,\text{d}y, \end{gather}

(2.17c)\begin{gather}{\bar{\hat{\mathcal{E}}}}={-}{\bar{\hat{\mathcal{W}}}} ={\rm d}\overline{{\hat\eta}_t^2} -{R}^{{-}1} [\overline{\hat\phi_y\hat\phi_{yy}} -3\overline{\hat\phi_{y}\hat\phi_{xx}}+U_y \bar{{\breve \phi}}_y]^{y=1}, \end{gather}

(2.17d)\begin{gather}{\bar{\hat{\mathcal{D}}}}={R}^{{-}1}\int_0^1 (4\overline{{\hat\phi}_{xy}^2}+ \overline{({\hat\phi}_{yy}-{\hat\phi}_{xx})^2}+2U_y \bar{{\breve\phi}}_{yy})\,{\rm d} y. \end{gather}

We particularly consider the averaged work done by the fluid on the wall over a period and wavelength, ${\bar {\hat {\mathcal {W}}}}=-{\bar {\hat {\mathcal {E}}}}$, which we compute below. This time-averaged energy budget (2.17) also serves as a useful check of our neutrally stable solutions and we confirm that the total error is $O(10^{-5})$ for all cases tested (Wang Reference Wang2019). Note also that this energy budget explicitly separates the work done by the fluid on the wall and the work done by the axial pressure gradient, whereas in Stewart et al. (Reference Stewart, Waters and Jensen2010b) these contributions were lumped into other terms.

2.7. Activation energy

As shown previously, the Reynolds–Orr energy budget calculation does not lead to a self-consistent prediction of the activation energy for a viscous perturbation mode. Instead, we follow the formulation of Landahl (Reference Landahl1962), who computed the activation energy of a neutrally stable mode as proportional to the difference in mechanical impedance arising from a small perturbation from the neutral point. As noted by Landahl (Reference Landahl1962), this approach does not necessarily hold for a viscous flow as it neglects some terms, but, as in that study, we assume that the sign of the activation energy is predicted correctly for a viscous flow. Indeed, we show below that the sign of the activation energy for the FISI is consistent with the classification framework of Benjamin (Reference Benjamin1960) based on the response to wall damping. Similarly, Cairns (Reference Cairns1979) defined activation energy for a neutral inviscid wave (using piecewise-constant velocity profiles with no shear or critical layers), arguing that since both boundary layer flow and plane Poiseuille flow have nearly neutral inviscid normal modes in the long-wave limit, their activation energy can be approximated using the inviscid equations (provided a small indentation is taken in the complex plane to pass around the critical point) and the sign of their activation energy predicts the destabilising effects of viscosity in the Orr–Sommerfeld equation. We follow the notation of Cairns (Reference Cairns1979) and express the activation energy in terms of the fluid pressure on the wall through the dispersion relation $D(\omega,k)$ (cf. (2.11d)); since ${D}=0$ for normal modes, it emerges that this difference in normal stress (impedance) must be Taylor expanded about the neutral stability point to obtain an expression for the activation energy in the form (Landahl Reference Landahl1962; Cairns Reference Cairns1979)

(2.18)\begin{equation} {\hat W} = {{\rm Re} }\left(\frac{1}{4}\omega_n \frac{\partial {{D}}}{\partial \omega}(k_n,\omega_n)\right) {\hat{\hat \eta}}^2 \equiv W {\hat{\hat \eta}}^2. \end{equation}

Based on the dispersion relation (2.11d), this formula suggests there are two contributions to the activation energy for a viscous mode at neutral stability, one arising due to viscous effects and the other due to wall elasticity (wall mass and pre-tension). However, we note that if ${{\phi }}$ is purely real then the viscous contribution to ${D}$ is purely complex, and so does not contribute to the activation energy. For an inviscid fluid (see Appendix A) there is another contribution which cancels here due to the no-slip boundary conditions. We use (2.18) below to compute the activation energy of neutrally stable FISI using a combination of numerical and asymptotic analysis. In particular, we apply (2.18) to compute explicit expressions for the activation energy in the limit of large wall mass (§ 3) and large wall damping (§ 4). Furthermore, to numerically compute the activation energy of a neutrally stable normal mode, we perturb the neutrally stable frequency by a real (small) amount $\pm {\rm \Delta} \omega$, and use a second-order centred finite difference stencil to approximate the derivative in (2.18) with error $O({\rm \Delta} \omega ^2)$, in the form

(2.19)\begin{equation} {W} \approx {\rm Re} \left(\frac{\omega_n}{4} \frac{{D}(\omega_n+{\rm \Delta} \omega,k_n) - {D}(\omega_n-{\rm \Delta} \omega,k_n)}{2 {\rm \Delta} \omega} \right). \end{equation}

This numerical interpretation is equivalent to computing the difference in mechanical impedance arising from a small perturbation from the neutral stability point, consistent with Landahl (Reference Landahl1962). Note this approach to construct the activation energy only works for normal modes which are surface based (flutter, divergence, TWF, SD) and does not work for the hydrodynamic modes (TS).

3.1. Numerical results

To assess the stability of viscous system in the absence of wall damping, in figure 2 we trace neutral stability curves as a function Reynolds number in terms of the critical wavenumber (figure 2a) and the corresponding wavespeed (figure 2b) for selected values of the wall mass. As expected, the system is unstable to class A TS modes within a tongue-shaped region for sufficiently large Reynolds numbers, analogous to those found in a rigid-walled channel (Drazin & Reid Reference Drazin and Reid1981). For small and large values of the wall mass the corresponding TS neutral stability curve is only modified very slightly from the case with zero wall mass (figure 2a). For example, the corresponding critical Reynolds number for instability is modified from $R_c\approx 7743$ for $m=0$ to $R_c \approx 7693$ for $m=1000$. However, we note that for a range of intermediate values of the wall mass the behaviour of the TS mode is more complicated, exhibiting a mode interaction with TWF which is explored in more detail in § 3.2.

Figure 2.

Neutral stability curves for Poiseuille flow as a function of the Reynolds number for selected values of the wall mass with fixed wall tension ($T=10$) and zero wall damping ($d=0$), illustrating: (a) critical wavenumber $k_n$; (b) critical wavespeed $c_n$. In each panel we consider $m=0$ (black), $m=1$ (red), $m=10$ (blue) and $m=100$ (magenta). The asterisks correspond to the leading-order prediction for neutrally stable TWF in the limit of large wall mass and low Reynolds number (3.2) while the open circles in (b) correspond to the leading-order wavespeed for large mass and large Reynolds number. The shaded regions in (a) indicate where the base state is asymptotically stable for $m=100$. The inset to (a) shows the neutrally stable wavenumber as a function of wall tension $T$ for three values of the Reynolds number ($R=1$, $R=10^3$ and $R=10^6$).

The base flow is also unstable to a single surface-based normal mode, which takes the form of class B TWF at neutral stability. The corresponding neutral stability curves each enclose a region for small wavenumbers and each exhibit two regimes (figure 2a,b). For small Reynolds numbers these neutral stability curves are independent of $R$, where the critical wavespeed reduces from twice the maximal flow speed ($c=3$ for $R=0$, where the system is notably still unstable to TWF) towards zero as the wall mass increases. Note that for our choice of non-dimensionalisation the limit $R\to 0$ corresponds to increasing fluid viscosity rather than vanishing flow velocity. Indeed, in a different approach where the wall parameters are scaled consistent with Squire's theorem (Rotenberry & Saffman Reference Rotenberry and Saffman1990), the TWF neutral stability curves no longer persist all the way to $R=0$ (Lucey & Carpenter Reference Lucey and Carpenter1995; Davies & Carpenter Reference Davies and Carpenter1997a).

For large Reynolds numbers and low wall mass, the neutrally stable TWF modes are long wavelength but the corresponding limiting wavespeed becomes less than the maximal flow speed ($c<3/2$) as the wall mass increases and the weak critical layer structure of Stewart et al. (Reference Stewart, Waters and Jensen2010b) is replaced by conventional critical layers (discussed in § 3.3). Conversely, for large values of the wall mass the critical wavenumber increases for large Reynolds numbers, encompassing modes of shorter wavelength (this short-wavelength mode is explored in detail in § 3.3). The transition between the long- and short-wavelength scalings occurs when $m\approx T$, reminiscent of the change in stability observed in potential flow with zero wall damping (Appendix A). However, in the viscous case the instability is class B and occurs as a result of a single normal mode becoming unstable rather than as a result of mode coalescence. Conversely, for $T>m$ the viscous system is long wavelength unstable while the potential flow system is stable. Hence, despite some qualitative similarity of the neutral stability curves, the flutter and TWF instabilities are distinct.

To further explore the behaviour of neutrally stable TWF as the wall becomes heavier, in figure 3 we trace the corresponding neutral stability curves as a function of wall mass for various Reynolds numbers, plotting the critical wavenumber (figure 3a) and the critical wavespeed (figure 3b) in the absence of wall damping. For low values of the wall mass the system is long wavelength unstable to TWF, where the neutral stability curves are independent of mass for small $m$ (figure 3a,b). Conversely, for large wall mass the critical wavenumber (figure 3a) and frequency (figure 3b) both increase with $m^{1/2}$ and the neutral stability curves become independent of Reynolds number; we exploit this asymptotic scaling in § 3.3. For the neutrally stable TWF mode the total work done by the fluid on the wall over a period is non-zero (figure 3c), in agreement with the mechanism of TWF described by Huang (Reference Huang1998). However, we note that this quantity changes sign when $m\approx T$, suggesting a change in mechanism (figure 3(c) inset). In particular, the total work done by the fluid normal stress on the wall is positive (${\bar {\hat {\mathcal {W}}}}=-{\bar {\hat {\mathcal {E}}}}$) for large Reynolds numbers and low wall mass. On the other hand, the activation energy of TWF is always positive for all values of wall mass, consistent with its classification as a class B instability (figure 3d). In summary, this figure shows that in the absence of wall damping the viscous system is unstable to class B TWF which always has positive activation energy but where the work done by the fluid on the wall can change sign as a function of wall mass.

Figure 3.

Neutral stability curves of TWF as a function of wall mass for various values of the Reynolds number with fixed pre-tension ($T=10$) and no wall damping ($d=0$): (a) critical wavenumber $k_n$; (b) critical wavespeed $c_n$; (c) the negative of the work done by the fluid stress on the wall $\bar {\hat {\mathcal {E}}}=-\bar {\hat {\mathcal {W}}}$; (d) the activation energy $W$. In each panel we consider $R=10$ (solid), $R=100$ (dashed) and $R=1000$ (dotted). The asterisks in (a,b,d) correspond to the leading-order prediction for neutrally stable TWF in the limit of large wall mass (3.2). The shaded region in (a) indicates where the base state is asymptotically stable for $R=1000$. The insets to (c) and (d) show the work done by the fluid on the wall and the activation energy on a linear scale for small values.

3.2. Mode interaction between TS and TWF

In § 3.1 we noted that as the wall mass increases the critical wavespeed of TWF reduces below the maximal flow speed and becomes independent of Reynolds number (figure 2). It emerges that for some parameter values the TWF mode and the TS mode interact (evident from the modification to the TS neutral stability curve in figure 2 for $m=100$). Such a mode interaction between class A and class B modes is energetically favourable.

To elucidate this mode interaction in more detail, figure 4 explores the neutral stability curves of TWF and TS modes in terms of the critical wavenumber (figure 4a) and wavespeed (figure 4b) as a function of the wall mass for fixed pre-tension and Reynolds number. Across most of the parameter space these neutral stability curves are distinct, but we note that for $50 \lesssim m \lesssim 120$ the neutral stability curve for TS is no longer independent of wall mass (this range is within the unstable zone of the primary TWF mode). Across this range the frequency of the TS and TWF modes become comparable (for a given wavenumber) and the resulting mode interaction breaks the TS neutral curve into two pieces. These two segments connect to the unstable zones for TS at large and small mass, respectively, but also exhibit a distinct loop at each end within which the base flow is asymptotically stable (shaded regions in figure 4). To examine these islands of stability, we plot the growth rate of the two most unstable normal modes (${\rm Im} (\omega )$) as a function of wavenumber for various values of the wall mass (figure 4c–e). For low mass (figure 4c) the TWF and TS modes remain distinct, and both are unstable across a narrow range of wavenumbers. The mode interaction is evident for $m=65$ (figure 4d), where the traces of the two modes deviate and almost touch one another, resulting in a narrow range of $k$ where both modes are stable within the larger TWF unstable zone, but also a substantial increase in the growth rate for larger wavenumbers. The mode interaction is similarly evident for $m=100$ (figure 4e), where the stable island is now relatively large compared with the outer TWF neutral curve (figure 4a). In summary, this figure demonstrates that TS and TWF modes can interact resulting in large stable islands within previously unstable regions.

Figure 4.

Mode interaction between TS and TWF, plotting the neutral stability curves for Poiseuille as a function of wall mass for fixed pre-tension ($T=10$) and Reynolds number ($R=10^6$) with no wall damping ($d=0$), illustrating: (a) the critical wavenumber $k_n$; (b) the critical wavespeed $c$. The shaded regions in (a) indicate where the base state is asymptotically stable. The asterisks in (a), (b) correspond to the leading-order prediction for neutrally stable TWF in the limit of large wall mass (3.2). To illustrate the mode interaction we plot traces of the perturbation growth rate ${\rm Im} (\omega )$ as a function of wavenumber for the two most unstable normal modes of the system for: (c) $m=10$; (d) $m=65$ and (e) $m=100$. The insets in (d) and (e) show a close-up of the mode interaction.

3.3. The limit of large wall mass

The viscous system with no wall damping exhibits an asymptotic scaling in the limit of large mass ($m\gg 1$) for fixed wall tension (figure 4), which we exploit by adopting the short-wavelength rescaling

(3.1a)\begin{equation} k = m^{1/2} k_0, \quad \omega =\omega_0+ m^{{-}1/2}\omega_1 + O(m^{{-}1}), \end{equation}

where $k_0$ and $\omega _0$ are to be determined at neutral stability. In this limit it emerges that the perturbation eigenfunction is zero across the interior of the channel and close to the lower rigid wall and exhibits a rapidly varying (viscous) solution over a narrow boundary layer adjacent to the compliant wall; a typical eigenfunction profile in this regime is shown in figure 5(a). Across this wall layer we expand the variables according to

(3.1b)\begin{equation} y =1- m^{{-}1/2} Y, \quad \phi(y) ={-} m^{{-}1/2} \varPhi_{0}(Y) + O(m^{{-}1}), \quad D=m D_0 + m^{1/2} D_1 + O(1). \end{equation}

In this limit, the inhomogeneous Orr–Sommerfeld system (2.11) takes the leading-order form

(3.1c)\begin{equation} \varPhi_{0,YYYY}-2k_0^2\varPhi_{0,YY}+k_0^4\varPhi_{0}=0, \quad (Y\ge 0), \end{equation}

subject to boundary conditions

(3.1d)\begin{equation} \varPhi_{0,Y}(0)={-}U_y(1)=6,\quad \varPhi_{0}(0)={-}\frac{\omega_0}{k_0}, \quad \varPhi_{0} \rightarrow 0 \ \text{as} \ Y\rightarrow \infty. \end{equation}

These equations can be solved explicitly in the form

(3.1e)\begin{equation} \varPhi_{0}={-}\left(\frac{\omega_0}{k_0}+ (\omega_0-6)Y\right)\exp({-}k_0Y), \quad (Y\ge 0). \end{equation}

For $k_0$ and $\omega _0$ as calculated below, this profile shows excellent agreement with the numerically computed perturbation streamfunction for large wall mass (figure 5a). As in the full system, the leading-order normal stress balance on the wall (2.11d) can be expressed as a dispersion relation

(3.1f)\begin{equation} D_0(\omega_0,k_0) \equiv \omega_0^2 - T k_0^2, \end{equation}

implying that $\omega _0=\sqrt {T} k_0$ for normal modes, and so the leading-order system is always neutrally stable. This is consistent with the potential flow formulation in Appendix A, where the inviscid class B mode is neutrally stable.

Figure 5.

Neutrally stable eigenfunctions for TWF for pre-tension $T=10$ and no wall damping, illustrating: (a) numerically computed eigenfunction profile ${\rm Re} (\phi (y))$ for $m=100$ and $R=10$ compared with the large $m$ approximation to the eigenfunction (3.1e) for the same wall mass; (b) numerically computed eigenfunction profile ${\rm Re} (\phi (y))$ for $m=10$, $R=10^6$. The dashed lines in (b) indicate the positions of the two critical layers. The inset to (b) shows the streamwise velocity profile ${\rm Re} (\phi _y(y))$ in the neighbourhood of the critical layer.

Expanding the normal stress balance (2.11d) to the following order ($O(m^{1/2})$), we obtain the next-order correction to the dispersion relation (which includes viscous effects) in the form

(3.1g)\begin{equation} D_1(\omega_0,\omega_1,k_0) \equiv \frac{1}{i k_0 R} (\varPhi_{0,YYY}(0)-k_0^2\varPhi_{0,Y}(0)) + 2 \omega_0\omega_1. \end{equation}

Normal modes require $D_1=0$, so that

(3.1h)\begin{equation} \omega_1 ={-}\frac{i k_0}{R}\frac{(\omega_0-6)}{\omega_0}, \end{equation}

independent of both the wall tension and the wall mass (in contrast to the inviscid flutter modes of Appendix A). Hence, neutrally stable TWF (with ${\rm Im} (\omega _1)=0$) requires $\omega _{0n}=6$ and $k_{0n}=6/\sqrt {T}$ in the limit of large wall mass, independent of the Reynolds number. We note particularly that the TWF can be shown to be unstable even for $R=0$ (unlike the TWF modes reported by Stewart et al. Reference Stewart, Waters and Jensen2010b) since the leading-order balance is between the tension and mass components of the flexible wall.

Expressing in terms of the original model parameters, we obtain the criticality conditions

(3.2)\begin{equation} k_n = 6\left(\frac{m}{T}\right)^{1/2}, \quad \omega_n = 6, \end{equation}

which show excellent agreement with the numerical predictions for large mass (asterisks in figures 2–4). The leading-order stability criterion involves only the viscous components of the flow and is independent of the fluid inertia. In this case, for $m\gg 1$ the instability is driven by the wall inertia alone, which dominates the leading-order activation energy (obtained by expanding (2.18)), in the form

(3.3)\begin{equation} W = \tfrac{1}{2}m\omega_0^2 + O(m^{1/2}), \end{equation}

the kinetic energy of wall deflections at leading order. Since ${{W}}>0$, this TWF instability is class B, as expected. This prediction shows almost perfect agreement with the numerical simulations for large wall mass (figure 3d).

At neutral stability the steady correction to the pressure gradient along the channel (2.13) takes the form

(3.4)\begin{equation} \bar{{\breve p}}_x = 72 R^{{-}1} + O(m^{{-}1/2}), \end{equation}

which remains $O(1)$ as $m\rightarrow \infty$. Hence, the breved terms are all $O(1)$ and do not feature in the leading-order ($O(m^{1/2}$)) components of the time-averaged energy budget (2.17); at neutral stability the time-averaged components take the simple form

(3.5)\begin{equation} \bar{\hat{\mathcal{D}}} =\bar{\hat{\mathcal{W}}}= \frac{288 k_{0n} m^{1/2}}{R} + O(1), \end{equation}

which is plotted in figure 3(c), showing excellent agreement with the numerical calculation.

Note that although asymptotically large values of $m$ are not feasible for conventional materials, the scalings show reasonable agreement with the numerical simulations even for $O(1)$ values of the wall mass (figure 3) and so this analysis provides insight into the mechanism of TWF in feasible parameter regimes. At leading order the oscillating wall and the fluid are perfectly in phase ($\phi$ is purely real when normalised by $\eta$), where the net work done by the fluid stress on the wall releases energy ($\bar {\hat {\mathcal {W}}}>0$) which is then dissipated by viscous effects across the boundary layer beneath the wall ($\bar {\hat {\mathcal {D}}}>0$). The neutrally stable TWF oscillation results in a slightly less negative (steady) pressure gradient along the channel compared with the baseline flow ($\overline {{\breve p}_x}>0$), but the accompanying decrease in energy is a higher-order effect i.e. the work done by Reynolds stresses does not feature at leading order in the energy balance. This observation contrasts with the TWF modes of Lebbal et al. (Reference Lebbal, Alizard and Pier2022), who found TWF modes driven by the working of Reynolds stresses, although the rate of working of fluid normal stress on the wall is also a (less significant) source of energy.

4.1. Numerical results

We consider the numerical results in two stages. First, we consider weak damping (§ 4.1.1) to examine how the class B TWF modes from § 3 are perturbed by an irreversible energy transfer to the wall. We then consider large wall damping (§ 4.1.2), where unstable SD modes destabilise for much lower Reynolds numbers and rearrange the structure of the neutral stability curves.

4.1.1. Weak wall damping

To examine the influence of weak wall damping on the TWF normal mode of § 3.1, figure 6 illustrates the TWF neutral stability curves as the wall damping is increased from low values, showing the critical wavenumber (figure 6a) and the critical frequency (figure 6b) as a function of the Reynolds number. TWF is stabilised by wall damping (as expected), where the critical wavenumber and frequency are both reduced for all values of the Reynolds number; the short-wavelength scaling of the neutral stability curve for large Reynolds numbers is overwhelmed by even small amounts of damping (see curve for $d=0.25$), with the critical wavenumber decreasing weakly with $R$ for finite damping.

Figure 6.

Neutrally stable curves of Poiseuille flow as a function of Reynolds number for various small values of the wall damping with fixed wall mass ($m=10$) and pre-tension ($T=10$) illustrating: (a) the critical wavenumber $k_n$ ; (b) the critical frequency $\omega _n$. In each panel we consider $d=0$ (solid line), $d=0.25$ (dashed line), $d=0.5$ (dashed line), $d=0.75$ (dash-dotted line), $d=0.9$ (dotted line) and $d=1$ (dash-dotted line). The shaded regions in (a) indicate where the base state is asymptotically stable for $d=0.5$. The insets indicate the neutral stability curves traced as a function of $R$ for the TS mode in terms of (a) critical wavenumber, (b) critical frequency. Normal modes of the weakly damped system including: (c) trace of the five most unstable eigenvalues in the complex frequency plane as a function of wavenumber for $R=10\ 000$ and $d=0.9$, where arrows indicate the direction of increasing wavenumber and red crosses indicate eigenvalues for $k=1$; (d) wider eigenvalue spectrum for $k=1$ and $R=10\ 000$, with an inset around the two most unstable centre modes.

As the wall damping increases a second neutral stability curve also becomes evident (see figure 6 for $d=0.5$, $d=0.75$, $d=0.9$), encompassing a stable island within the TWF unstable region. As the wavenumber increases, the locus of the surface-based mode in the complex frequency plane is deviated toward the locus of another hydrodynamic mode (figure 6c), causing it to cross into the lower half-plane and form a stable zone. We identify this hydrodynamic mode by continuing to larger wavenumbers, showing that it takes the form of a ‘centre’ mode, located on the P branch of the eigenvalue spectrum with wavespeed close to the maximal flow speed (Schmid & Henningson Reference Schmid and Henningson2001); see the corresponding eigenvalue spectrum in figure 6(d). Such centre modes have recently been identified as destabilising in channel flows of viscoelastic fluids (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021), but this mode always remains stable in the present case. This stable island enlarges as the wall damping increases, where the encompassing neutral stability curve interacts with the outer TWF neutral stability curve, deviating toward one another (e.g. see neutral stability curves for $d=0.9$ in figure 6a,b). For sufficiently large wall damping this stable loop interacts with the TWF neutral curve, merging and separating in a different configuration (see curves for $d=1$ in figure 6a,b); this neutral curve structure is discussed in greater detail in § 4.1.2. Note that the neutral stability curves for TS (insets to figure 6a,b) are not significantly influenced by wall damping. In summary, this figure demonstrates that neutrally stable TWF is stabilised by wall damping, reducing the critical wavenumber but also facilitating a mode interaction which creates a stable island within the TWF unstable region. Furthermore, figure 6 shows how the three neutrally stable points for a given Reynolds number (for sufficiently large wall damping) all arise from instability of the same normal mode.

4.1.2. Large wall damping

To examine the stability of the system for values of the wall damping beyond the merger noted previously, figure 7 illustrates the neutral stability curves as a function of Reynolds number in terms of the critical wavenumber (figure 7a) and the critical frequency (figure 7b). The neutral stability curves are rearranged compared with figure 6: the low Reynolds number portion of the TWF neutral curve has merged with the lower branch of the stable loop, forming a single neutral stability curve for TWF, while the large-Reynolds-number portion of the TWF neutral stability curve has merged with the upper branch of the stable loop forming a two branch loop for class A SD. Hence, similar to Stewart et al. (Reference Stewart, Waters and Jensen2010b), SD is only unstable above a critical Reynolds number and the neutral stability curve exhibits a two branch structure; we term the upper (lower) branch as that with greater (lower) critical wavenumber. For non-zero wall damping the time-averaged work done by the fluid normal stress on the wall ($\bar {\hat {\mathcal {W}}}$) is always positive for all normal modes (figure 7c). We note that along the upper branch of the SD neutral stability curve the activation energy is weakly negative (figure 7d), approaching zero from below as the Reynolds number increases. This suggests that the upper branch SD for large damping is nearly class C, reminiscent of the class C divergence mode that becomes unstable in damped potential flow (Appendix A). However, we show in § 4.2 that this upper-branch neutral stability curve approaches zero wavenumber as the Reynolds number increases, and so the critical wavenumber does not asymptote to $k_d$ as $R \rightarrow \infty$. We also show below that this upper-branch SD mode exhibits a very different energy balance to the inviscid divergence mode, involving viscous effects at leading order, and so we cannot conclude that these two modes are aligned.

Figure 7.

Neutral stability curves of Poiseuille flow as a function of Reynolds number for various large values of the wall damping with fixed wall mass ($m=10$) and pre-tension ($T=10$), illustrating: (a) the critical wavenumber $k_n$; (b) the critical frequency $\omega _n$; (c) the negative of the work done by the fluid stress on the wall $\bar {\hat {\mathcal {E}}}=-\bar {\hat {\mathcal {W}}}$; (d) the activation energy $W$. In each panel we consider $d=10$ (black), $d=100$ (red) and $d=1000$ (blue). The shaded region in (a) indicates where the base state is asymptotically stable for $d=10$. The symbols in (a), (b) and (d) correspond to the leading-order prediction of the large $d$ approximation for TWF (asterisks, (4.13a,b) with $+$), SD(l) (filled triangles, (4.13a,b) with $-$) and SD(u) (open inverted triangles, (4.38a,b)). The insets in (a) and (b) show the neutral stability curves as a function of Reynolds number for the large $d$ theory (filled circles, (4.37a,b)). The inset to (c) shows the ratio of the work done by fluid normal stress on the wall ($\bar {\hat {\mathcal {W}}}$) to the energy dissipation in the bulk fluid ($\bar {\hat {\mathcal {D}}}$) computed numerically for SD as a function of Reynolds number, alongside asymptotic predictions in the limit of large damping for both lower-branch SD (inverted triangles, computed from (4.16)) and upper-branch SD (triangles, computed numerically from (4.19)). The inset in (d) shows a close-up around zero activation energy.

In order to assess the properties of these FISI, figure 8 plots neutral stability curves of the system as a function of wall damping. As expected, for each Reynolds number the system exhibits three neutral stability curves for sufficiently large wall damping; the neutral stability curves for TWF and lower branch SD are parallel and exhibit similar critical wavenumber (figure 8a) and frequency (figure 8b), where both scale with $d^{-1/2}$ as $d\rightarrow \infty$. Conversely, upper-branch SD exhibits critical wavenumber independent of wall damping (figure 8a) and the critical frequency scales with $d^{-1}$ as $d \rightarrow \infty$ (figure 8b). We utilise these scalings in § 4.2 to construct rational asymptotic descriptions of TWF and SD in the limit of large wall damping.

Figure 8.

Neutral stability curves of Poiseuille flow as a function of the wall damping for various Reynolds numbers with fixed wall mass ($m=10$) and pre-tension ($T=10$), illustrating: (a) the critical wavenumber $k_n$; (b) the critical frequency $\omega _n$. In each panel we consider $R=10$ (black), $R=100$ (red) and $R=1000$ (blue). The shaded region in (a) indicates where the base state is asymptotically stable for $R=10$. The symbols correspond to the leading-order prediction of the large $d$ approximation for TWF (asterisks, (4.13a,b) with $+$), SD(l) (filled triangles, (4.13a,b) with $-$) and SD(u) (open inverted triangles, (4.38a,b)).

To elucidate the spatial structure of the three neutrally stable FISI, figure 9 illustrates the streamfunction profiles of TWF (figure 9a), lower-branch SD (figure 9b) and upper-branch SD (figure 9c) for the same Reynolds number and choice of wall properties. The spatial profiles of TWF and lower-branch SD are similar, exhibiting a non-trivial asymmetric profile across the entire domain. Conversely, upper-branch SD exhibits rapid variations across narrow boundary layers adjacent to the channel walls, indicating a substantial difference in mechanism.

Figure 9.

Neutrally stable eigenfunctions $\phi$ for (a) TWF; (b) SD(l); (c) SD(u). The asterisks in (a,b) correspond to the leading-order asymptotic prediction for large wall damping (4.6). The asterisks in (c) correspond to the inviscid core solution (4.22), while the dashed lines correspond to the two boundary layer solutions (4.26) and (4.30). Here, $T=10$, $m=10$ and $R=1000$.

4.2. The limit of large wall damping

The scalings identified in § 4.1 facilitate a rational asymptotic description of neutrally stable TWF and SD modes in the limit of large wall damping, building on the ad hoc description of lower-branch SD provided by Stewart et al. (Reference Stewart, Waters and Jensen2010b) (based on a flow profile assumption). In particular, in the limit of large wall damping we describe neutrally stable lower-branch SD and TWF (§ 4.2.1) and neutrally stable upper-branch SD (§ 4.2.2).

4.2.1. Neutrally stable TWF and lower-branch SD in the limit of large wall damping

Following the scalings identified in § 4.1 we return to the inhomogeneous system (2.11) and treat the wall damping as a large parameter ($d\gg 1$). We adopt a long-wavelength, low-frequency rescaling in the form

(4.1)\begin{equation} k = d^{{-}1/2} k_0, \quad \omega = d^{{-}1/2} \omega_0 + d^{{-}1} \omega_1 + O(d^{{-}3/2}), \end{equation}

where $k_0$ and $\omega _0$ are to be determined at neutral stability. We expand the perturbation streamfunction and the dispersion relation $D$ according to

(4.2)\begin{equation} \phi = \phi_0 + d^{{-}1/2} \phi_1 + O(d^{{-}1}), \quad D=d^{1/2} D_0 + D_1 + O(d^{{-}1/2}). \end{equation}

The inhomogeneous Orr–Sommerfeld system (2.11) becomes at leading order

(4.3)\begin{equation} \phi_{0,yyyy}=0, \quad (0\le y \le 1), \end{equation}

subject to boundary conditions

(4.4)\begin{gather} \phi_0=0, \quad \phi_{0,y}=0 \quad (y=0), \end{gather}

(4.5)\begin{gather}\phi_0=\frac{\omega_0}{k_0}, \quad \phi_{0,y}={-}U'(1) = 6, \quad (y=1). \end{gather}

We solve for $\phi _0$ explicitly in the form

(4.6)\begin{equation} \phi_0(y) = \frac{1}{k_0}y^2(\omega_0 (3 - 2y) + 6 k_0 ({-}1 + y)), \quad (0\le y \le 1). \end{equation}

Given $k_0$ and $\omega _0$ as calculated below, this profile is plotted in figure 9(a) (figure 9(b)) against numerical predictions of the eigenfunction for TWF (lower-branch SD) showing excellent agreement. In both cases the leading-order streamfunction is purely real for neutrally stable solutions.

Substituting these expansions into the normal stress balance (2.11d) and substituting for $\phi _0$ we obtain the leading-order dispersion relation

(4.7)\begin{equation} D_0(\omega_0,k_0) \equiv \frac{1}{ik_0 R}\phi_{0,yyy}(1) + i\omega_0 = \frac{12}{ik_0 R}\left(\frac{3\omega_0}{k_0}-1\right) + i\omega_0, \end{equation}

where the first term corresponds to viscous effects and the second term is due to wall damping. For normal modes $D_0(\omega _0,k_0)=0$, so we have

(4.8)\begin{equation} \omega_0= \frac{36 k_0}{12 + k_0^2 R}, \end{equation}

where $\omega _0$ is purely real for $k_0$ real.

To fix the corresponding critical wavenumber of the system we expand the governing equations to the following order in the form:

(4.9a)\begin{equation} \frac{1}{ik_0 R}{\phi}_{1,yyyy} = \left(U - \frac{\omega_0}{k_0}\right) {\phi}_{0,yy} -U_{yy}{\phi}_0, \end{equation}

subject to the inhomogeneous boundary conditions

(4.9b)\begin{gather} {\phi}_1(0)=0, \quad {\phi_{1,y}}(0)=0, \end{gather}

(4.9c)\begin{gather}{\phi}_1(1)=\frac{\omega_1}{k_0}, \quad {\phi_{1,y}}(1)=0. \end{gather}

We integrate (4.9) four times and apply the four boundary conditions to determine $\phi _1$ explicitly as a polynomial in $y$. The normal stress condition (2.11d) at $O(1)$ then implies the first-order correction to the dispersion relation, in the form

(4.10)\begin{equation} D_1(\omega_0,\omega_1,k_1) \equiv \frac{1}{i k_0 R}\phi_{1,yyy}(1) + i\omega_1. \end{equation}

Setting $D_1=0$ and substituting for $\phi _1(y)$, we obtain the correction to the frequency in the form

(4.11)\begin{equation} \omega_1 = \frac{i \omega_0 R}{210 k_0}(9k_0^2 - 17\omega_0 k_0 + 7 \omega_0^2 ). \end{equation}

Looking for neutrally stable solutions (${\rm Im} (\omega _1)=0$), we obtain a second constraint linking $k_0$ and $\omega _0$ in the form

(4.12)\begin{equation} 9k_0^2 - 17\omega_0 k_0 + 7 \omega_0^2=0. \end{equation}

Solving the constraints (4.8) and (4.12) simultaneously, we obtain a pair of physical solutions which can be rewritten in terms of the unscaled dimensionless variables and parameters in the form

(4.13a,b)\begin{equation} k_{n\pm} = \left(\frac{2(11 \mp \sqrt{37})}{dR} \right)^{1/2}, \quad \omega_{n\pm} = \frac{2}{7}\left(\frac{291\pm 6 \sqrt{37}}{dR} \right)^{1/2}. \end{equation}

These solutions are almost indistinguishable from the numerically computed neutral stability curves as $d$ becomes large (figures 7a,b and 8a,b), where $k_{n+}$ corresponds to TWF and $k_{n-}$ corresponds to lower-branch SD.

The dispersion relation (2.11d) expanded to first order ($D\approx d^{1/2}D_0+D_1$) can be differentiated to calculate the leading-order activation energy (2.18). Applying the neutral stability condition (4.8), the corresponding activation energy can be calculated as

(4.14)\begin{equation} {{W}}_\pm = {\frac{3 (37 \pm 17 \sqrt{37})}{980}} + O(d^{{-}1/2}), \end{equation}

independent of Reynolds number. So neutrally stable TWF (SD), the $+(-)$ branch of solutions, has positive (negative) activation energy, as expected. Again these show excellent agreement with numerically computed activation energy along the two neutrally stable branches (figure 7d).

This asymptotic analysis elucidates the steady correction to driving pressure gradient (2.13), which can be calculated as

(4.15)\begin{equation} \bar{\breve p}_x ={-}\frac{36(25 \pm \sqrt{37})}{7R} + O(d^{{-}1/2}). \end{equation}

Hence, unlike TWF modes described in the limit of large mass (§ 3.3), the oscillating flow serves to make the steady pressure gradient along the channel more negative at leading order.

Furthermore, we calculate the leading-order contributions to the terms in the time-averaged energy budget (2.17), where we obtain

(4.16)\begin{equation} (\bar{\hat{\mathcal{D}}},\bar{\hat{\mathcal{W}}},\bar{\hat{\mathcal{P}}})=\frac{(25\pm\sqrt{37})}{49R}(804,552,252) +O(d^{{-}1/2}). \end{equation}

Hence, approximately 2/3 of the energy consumed by dissipation across the bulk ($\bar {\hat {\mathcal {D}}}$) is provided by the working of fluid stresses on the wall ($\bar {\hat {\mathcal {W}}}$), in strong agreement with the numerical simulations (inset to figure 7c), with the remaining amount released by a decrease in the steady (negative) pressure gradient along the channel ($\bar {\hat {\mathcal {P}}}$). However, in this case we also note that the contributions involving breved variables to the leading-order terms in the time-averaged energy budget (2.17) combine perfectly to zero at this order, and so there is an additional balance between the terms involving hatted variables in $\bar {\hat {\mathcal {W}}}$ and $\bar {\hat {\mathcal {D}}}$, suggesting a similar mechanism to the TWF mode described in § 3.3 where the correction to the mean flow does not play a role; this balance is not true in general and arises because the leading-order eigenfunction is purely real when normalised on the wall deflection, so the wall and flow are perfectly in phase at neutral stability.

Note that, although asymptotically large values of wall damping are not feasible for conventional materials, the scalings show reasonable agreement with the numerical simulations even for $O(1)$ values (figure 8) and so this analysis provides insight into the mechanism of these FISI in feasible parameter regimes. At leading order the oscillating wall and the fluid are perfectly in phase ($\phi$ is purely real when normalised by $\eta$), where both the net work done by the fluid stress on the wall ($\bar {\hat {\mathcal {W}}}>0$) and the work done by the steady pressure gradient ($\bar {\hat {\mathcal {P}}}>0$) release energy (in a ratio of approximately 2 : 1), which is then dissipated by viscous effects across the channel ($\bar {\hat {\mathcal {D}}}>0$). A similar mechanism of SD driven by the working of fluid normal stress on the wall was recently noted by Lebbal et al. (Reference Lebbal, Alizard and Pier2022). However, their work found that TWF modes were primarily driven by the rate of working of Reynolds stresses, in contrast to the modes analysed here, although the working of fluid normal stress on the wall was also shown to be mildly destabilising.

4.2.2. Neutrally stable upper-branch SD in the limit of large wall damping

Following the scalings identified in figure 8 we return to the full model (2.11) and assume the wavenumber remains $O(1)$ with respect to the wall damping and rescale according to

(4.17)\begin{equation} k=k_0, \quad \omega = d^{{-}1} \omega_0 + O(d^{{-}2}), \end{equation}

where we calculate the criticality conditions in terms of the wavenumber $k_0$ and leading-order frequency $\omega _0$. Similarly, we expand the streamfunction according to

(4.18)\begin{equation} \phi = \phi_0 + O(d^{{-}1}). \end{equation}

The inhomogeneous linear stability problem (2.11) becomes at leading order

(4.19a)\begin{equation} U (\phi_{0,yy} - k_0^2 \phi_0) - U_{yy}\phi_0 = \frac{1}{i k_0 R}(\phi_{0,yyyy} - 2k_0^2 \phi_{0,yy}+ k_0^4 \phi_0), \end{equation}

subject to boundary conditions

(4.19b)\begin{gather} \phi_0=0, \quad \phi_{0,y}=0 \quad (y=0), \end{gather}

(4.19c)\begin{gather}\phi_0=0, \quad \phi_{0,y}={-}U_y(1) = 6, \quad (y=1); \end{gather}

note that the leading-order kinematic boundary condition does not involve the wall profile in this limit. This system for $\phi _0$ is notably independent of the leading-order frequency $\omega _0$. Substituting these expansions into the normal stress balance (2.11d) leads to the dispersion relation

(4.19d)\begin{equation} D(\omega_0,k_0) \equiv \frac{1}{ik_0 R}({\phi}_{0,yyy}(1) - k_0^2 \phi_{0,y}(1) ) - T k_0^2 + i \omega_0 + O(d^{{-}1}). \end{equation}

In this case, the wall mass does not feature in the leading-order problem, but, unlike lower-branch SD modes (§ 4.2.1), the restoring force due to wall tension features prominently. To the best of our knowledge it is not possible to solve this system analytically, but numerical solution is straightforward using a spectral method analogous to that described in § 2.4. It emerges that numerical solution of the leading-order system (4.19) is sufficient to isolate the critical conditions for instability; these neutral stability curves are shown in figure 7, plotting the approximated critical wavenumber (inset to figure 7a) and critical frequency (inset to figure 7b) as a function of the Reynolds number, showing excellent agreement with the numerically computed neutral stability curves of upper-branch SD.

Given the leading-order dispersion relation (4.19d) we find that the leading-order activation energy for neutrally stable upper branch SD is $O(d^{-1})$ and so approaches zero as the wall damping becomes large, consistent with the numerical simulations (figure 8d). Hence, this upper-branch SD mode approaches class C as the wall damping becomes large, reminiscent of the inviscid divergence mode (Appendix A).

This approach does not provide analytical insight into the partition of energy, but we note that numerical calculation of the corresponding energy terms in (2.17), plotted as symbols in figure 7(c), show excellent agreement with the full computations. These numerical computations indicate that the mechanism of upper-branch SD is broadly similar to lower-branch SD, with energy released by the working of fluid stress on the wall ($\bar {\hat {\mathcal {W}}}>0$) and the working of the steady pressure gradient along the channel ($\bar {\hat {\mathcal {P}}}$), which is then dissipated by viscous effects across the bulk ($\bar {\hat {\mathcal {D}}}$). In the limit of large damping the ratio of work done by fluid normal stress on the wall compared with the work done by dissipation in the fluid $\bar {\hat {\mathcal {W}}}/\bar {\hat {\mathcal {D}}} \approx 0.73$ as $R \rightarrow \infty$ (almost perfectly independent of $T$); this value is close to, but slightly larger than, that predicted for lower-branch SD ($\bar {\hat {\mathcal {W}}}/\bar {\hat {\mathcal {D}}} = 46/67 \approx 0.69$, § 4.2.1), but the balance is broadly similar. Interestingly, unlike lower-branch SD, the contributions to these energies which arise from the (breved) terms do not all cancel internally, and so contribute to the overall partition (the leading-order eigenfunction is complex, so the wall and streamfunction profiles are not in phase, and the work done by Reynolds stresses is non-negligible).

To form a more intuitive description of the flow we consider the leading-order system in the limit of large damping (4.19) in the additional limit of large Reynolds number. It emerges that the eigenfunction $\phi _0$ can be separated into an inviscid core with viscous effects confined to a narrow boundary (critical) layers adjacent to the walls. Based on the scalings evident in figure 7(a,b), we consider long-wavelength, low-frequency instabilities of the form

(4.20)\begin{equation} k_0 = R^{{-}1/7} k_{00}, \quad \omega_0 = R^{{-}2/7} \omega_{00} + O(R^{{-}4/7}). \end{equation}

Across the inviscid core, the eigenfunction $\phi _0$ is expanded according to

(4.21)\begin{equation} \phi_0 = \phi_{00} + R^{{-}2/7} \phi_{01} + O(R^{{-}4/7}), \end{equation}

where the governing equations (4.19) become at leading order in $R$

(4.22a)\begin{equation} U\phi_{00,yy} - U_{yy}\phi_{00} = 0, \end{equation}

subject to reduced boundary conditions

(4.22b)\begin{equation} \phi_{00}(0)=0, \quad \phi_{00}(1)=0. \end{equation}

This system admits the simple solution $\phi _{00} = A U$, where $A$ is an undetermined constant. Note that this inhomogeneous system has already been normalised with respect to the wall deflection, and so $A$ must be calculated as part of the solution. Plotting this parabolic profile (for $A$ as computed below) shows reasonable agreement with the numerically computed eigenfunction (figure 9c). In this case the outer limit of the axial velocity component takes the form

(4.22c)\begin{gather} \phi_{00,y} \rightarrow A U_y(0) = 6A \ \text{as} \ y\rightarrow 0, \end{gather}

(4.22d)\begin{gather}\phi_{00,y} \rightarrow A U_y(1) ={-}6A \ \text{as} \ y\rightarrow 1. \end{gather}

For large Reynolds numbers the flow exhibits narrow boundary layers adjacent to the walls of width $R^{-2/7}$. Adjacent to the rigid wall (denoted with the superscript $^{(r)}$) we rescale

(4.23)\begin{equation} y=R^{{-}2/7} Y, \quad \phi_0 = R^{{-}2/7} \varPhi^{(r)}_{00}(Y) + O(R^{{-}4/7}). \end{equation}

In this case, the governing equations (4.19) take the leading-order form across the boundary layer

(4.24a)\begin{equation} 6 Y k_{00} \varPhi^{(r)}_{00,YY} + i \varPhi^{(r)}_{00,YYYY} = 0, \quad (Y \ge 0), \end{equation}

subject to boundary conditions

(4.24b)\begin{equation} \varPhi^{(r)}_{00}(0)=0, \quad \varPhi^{(r)}_{00,Y}(0)=0, \quad \varPhi^{(r)}_{00,Y} \rightarrow 6A \ \text{as} \ Y \rightarrow \infty. \end{equation}

This is Airy's equation in $\varPhi ^{(r)}_{00,YY}$, with solution

(4.25)\begin{equation} \varPhi^{(r)}_{00,YY} = A_0 \text{Ai}((6ik_{00})^{1/3} Y ) + B_0 \text{Bi}((6ik_{00})^{1/3} Y ), \end{equation}

where $\text {Ai}$ and $\text {Bi}$ are Airy functions of the first and second kinds, respectively, and $A_0$ and $B_0$ are complex constants; we immediately set $B_0=0$ to ensure the solution is bounded at $Y=0$. Integrating this expression once, and applying the no-slip and far-field conditions (4.24b) as well as known results for the integral of Airy functions (Abramowitz & Stegun Reference Abramowitz and Stegun1965), we deduce

(4.26)\begin{equation} \varPhi^{(r)}_{00,Y}(Y) = 18 A (6ik_{00})^{1/3} \int_0^Y \text{Ai}((6ik_{00})^{1/3} \zeta)\,\text{d}\zeta. \end{equation}

For $k_{00}$ and $A$ calculated below, this profile is plotted in figure 9(c), showing good agreement with the numerically computed streamfunction profile close to the rigid wall. This expression can be integrated once more to obtain the streamfunction explicitly, although the final expression is cumbersome and not reproduced here. In the limit as $Y\rightarrow \infty$ we obtain the limiting form

(4.27a)\begin{equation} \varPhi^{(r)}_{00} \rightarrow \frac{18 A}{(6ik_{00})^{1/3}}\left( \frac{\varGamma(\frac{1}{3})(6ik_{00})^{1/3}Y}{9 \varGamma(\frac{4}{3})} - \frac{1}{3^{1/3} \varGamma(\frac{1}{3})} + O(Y^{{-}1})\right), \end{equation}

where $\varGamma ({\cdot })$ is the gamma function; this ejection of fluid from the boundary layer drives flow in the inviscid core at first-order. The leading-order fluid pressure is constant across the boundary layer on the rigid wall and takes the form

(4.27b)\begin{equation} P_{00}^{(r)} = \frac{1}{i k_{00}}\varPhi^{(r)}_{00,YYY}(0) ={-}\frac{18\times 6 A}{3^{1/3} (6i k_{00})^{1/3} \varGamma(\frac{1}{3})}. \end{equation}

Similarly, adjacent to the compliant wall (denoted with the superscript $^{(c)}$) we rescale the variables according to

(4.28)\begin{equation} y=1- R^{{-}2/7} Y, \quad \phi_0 = R^{{-}2/7} \varPhi^{(c)}_{00}(Y) + O(R^{{-}4/7}). \end{equation}

The governing equations (4.19a) across the boundary layer again reduce to Airy's equation (4.24a), subject to boundary conditions

(4.29)\begin{equation} \varPhi^{(c)}_{00}(0)=0, \quad \varPhi^{(c)}_{00,Y}(0)=U_y(1) ={-}6, \quad \varPhi^{(c)}_{00,Y} \rightarrow -A U_y(1) = 6A \ \text{as} \ Y \rightarrow \infty. \end{equation}

In a similar manner to the other boundary layer, we deduce

(4.30)\begin{equation} \varPhi^{(c)}_{00,Y}(Y) = U_y(1) + 18 (1+ A ) (6ik_{00})^{1/3} \int_0^Y \text{Ai}((6ik_{00})^{1/3} \zeta)\,\text{d}\zeta. \end{equation}

For $k_{00}$ and $A$ calculated below, this profile is also shown in figure 9(c), again showing close agreement to the numerically computed profile close to the flexible wall. As before, this expression can be integrated once more to obtain the limiting form of the streamfunction as $Y\rightarrow \infty$ in the form

(4.31a)\begin{equation} \varPhi^{(c)}_{00} \rightarrow -6Y + \frac{18 (A+1)}{(6ik_{00})^{1/3}}\left( \frac{\varGamma(\frac{1}{3}) (6ik_{00})^{1/3}Y}{9 \varGamma(\frac{4}{3})} - \frac{1}{3^{1/3} \varGamma(\frac{1}{3})} + O(Y^{{-}1})\right), \end{equation}

which, similar to the boundary layer on the rigid wall, drives flow in the inviscid core at first order. The leading-order fluid pressure is constant across the boundary layer on the compliant wall and takes the form

(4.31b)\begin{equation} P_{00}^{(c)} ={-}\frac{1}{i k_0}\varPhi^{(c)}_{00,YYY}(0) = \frac{18\times 6 (A+1)}{3^{1/3} (6 i k_{00})^{1/3} \varGamma(\frac{1}{3})}. \end{equation}

The upper and lower boundary layer profiles can be solved for a given $A$. To determine the constant $A$ we expand the inviscid core flow to the following order in $R^{-2/7}$ to obtain the governing equations

(4.32a)\begin{equation} U\phi_{01,yy} - U_{yy}\phi_{01} = (U\phi_{01,y} - U_y \phi_{01})_y= A k_{00}^2 U^2, \end{equation}

subject to boundary conditions modified by the flow ejected from the boundary layers (4.27a), (4.31a), in the form

(4.32b)\begin{equation} \phi_{01}(0)= \frac{18 A}{3^{1/3}(6ik_{00})^{1/3} \varGamma(\frac{1}{3})}, \quad \phi_{01}(1)= \frac{18 (A+1)}{3^{1/3}(6ik_{00})^{1/3}\varGamma(\frac{1}{3})}. \end{equation}

This system can be solved directly using the method of variation of parameters. However, we note that the $x$-momentum equation (2.6a) expanded to $O(R^{-2/7})$ across the inviscid core can be expressed as

(4.33)\begin{equation} p_{01} ={-}(U\phi_{01,y} - U_y \phi_{01}), \end{equation}

so that (4.32a) can be written in terms of the fluid pressure, in the form

(4.34a)\begin{equation} p_{01,y} ={-}A k_{00}^2 U^2. \end{equation}

Furthermore, using (4.33), the boundary conditions on $\phi _{01}$ (4.32b) can be re-written in terms of the fluid pressure in the lower (4.27b) and upper boundary layers (4.31b), such that

(4.34b)\begin{equation} p_{01}(0) = P_{00}^{(r)} ={-}\frac{18\times 6 A}{3^{1/3} (6i k_{00})^{1/3} \varGamma(\frac{1}{3})}, \quad p_{01}(1) = P_{00}^{(c)} = \frac{18\times 6 (A+1)}{3^{1/3} (6 i k_{00})^{1/3} \varGamma(\frac{1}{3})}. \end{equation}

The two boundary layers are interacting across the core flow, in a similar manner to the large-Reynolds-number approximation to lower-branch TS instabilities in a rigid channel (Bogdanova & Ryzhov Reference Bogdanova and Ryzhov1983; Pedley Reference Pedley2000; Stewart et al. Reference Stewart, Waters and Jensen2010b). Integrating (4.34a) once and applying both boundary conditions (4.34b), determines the amplitude of the leading-order flow in the inviscid core, in the form

(4.35)\begin{equation} A=\left({-}2 + \frac{6}{5}k_{00}^2\frac{3^{1/3} (6ik_{00})^{1/3} \varGamma(\frac{1}{3})}{108} \right)^{{-}1}. \end{equation}

We note that for $k_{00}\ll 1$ in (4.34a), where the cross-stream pressure gradient is weak, $P_{00}^{(r)}=P_{00}^{(c)}$ and this constraint reduces to $A=-\frac {1}{2}$.

The compliant wall boundary layer is coupled to the leading-order dispersion relation

(4.36)\begin{equation} D_1(\omega_{01},k_{00})\equiv P_{00}^{(c)} - T k_{00}^2 + i\omega_{00} = \frac{108 (A+1)}{3^{1/3} (6i k_{00})^{1/3} \varGamma(\frac{1}{3})} - T k_{00}^2 + i\omega_{00}. \end{equation}

Neutrally stable modes (${\rm Im} (\omega _{00})=0$) can be identified by setting $D_1=0$ coupled to the criterion for $A$ (4.35), resulting in a polynomial dispersion relation which can be solved for $k_{00n}$ and $\omega _{00n}$ by separating into real and imaginary parts; however, this final step must be taken numerically. In the simplified case where the cross-stream pressure gradient is weak we can approximate $A\approx -\frac {1}{2}$, these criticality conditions reduce to the explicit form

(4.37a,b)\begin{equation} k_{00n} = \left(\frac{3^{17/6}}{2^{1/3} T \varGamma(\frac{1}{3})} \right)^{3/7} , \quad \omega_{00n}=\frac{1}{3^{1/2}}T k_{00n}^2. \end{equation}

Furthermore, these criticality conditions can be expressed in terms of the unscaled dimensionless variables and parameters in the form

(4.38a,b)\begin{equation} k_{n} = \left(\frac{3^{17/6}}{2^{1/3} \varGamma(\frac{1}{3})} \right)^{3/7} T^{{-}3/7} R^{{-}1/7}, \quad \omega_{n}=\left(\frac{3^{9/4}}{2^{1/3} \varGamma(\frac{1}{3})} \right)^{6/7} T^{1/7} d^{{-}1} R^{{-}2/7}. \end{equation}

These expressions demonstrate excellent agreement with the numerically computed neutral stability curves shown in figures 7 and 8.

Note that this analytical calculation could be extended to explicitly compute the partition of energies for the neutrally stable modes (similar to lower-branch TWF and SD), where each term is formally $O(R^{-1})$. However, this calculation requires the first correction to the streamfunction across the viscous boundary layers, and so is not pursued further here.

4.3. Mode interaction with TS for large wall mass

In § 3.2 we elucidated a mode interaction between the undamped TWF and TS modes for increasing wall mass, exhibiting a new stable island within the unstable TWF zone. To explore this mode interaction in the presence of wall damping, in figure 10 we plot the critical wavenumber (figure 10a) and frequency (figure 10b) of the system as a function of Reynolds number for several values of the wall mass holding all other parameters fixed. For low wall mass (e.g. $m=10$ in figure 10), the wall damping is sufficiently large so that the SD neutral curve has disconnected from the TWF neutral curve (cf. § 4.2). As the wall mass increases the interaction with the centre mode, which leads to the separate TWF and SD neutral stability curves (§ 4.1), becomes weaker, meaning that the locus of the surface-based mode as a function of wavenumber no longer transiently crosses into the lower half-plane. Instead, the system exhibits a single neutral stability curve and the SD and TWF neutral curves merge together ($m=100$ in figure 10). As the wall mass increases the TS and TWF modes interact and eventually merge into a configuration reminiscent of the case with low wall damping, where the system exhibits a stable loop within the outer unstable region (e.g. $m=200$). For even larger values of the wall mass the neutral stability curve self-intersects, creating a second stable island in the parameter space (e.g. $m=1000$). Hence, although increasing wall mass is destabilising to TWF, mode interaction leads to new stable islands within the unstable region. Figure 10 also exhibits further evidence of mode interaction with a visible modification to the TS neutral stability curve for large wall mass (e.g. $m\geqslant 100$). However, given that this interaction occurs for very large Reynolds numbers we do not consider this further.

Figure 10.

Neutral stability curves for Poiseuille flow as a function of wall mass for various Reynolds numbers with fixed wall damping ($d=10$) and pre-tension ($T=10$), illustrating: (a) the critical wavenumber $k_n$; (b) critical frequency $\omega _n$. Here, we consider $m=10$ (solid line), $m=100$ (dashed line), $m=200$ (dash-dotted line), $m=1000$ (dotted line) and $m=10\ 000$ (dash-dotted line). The shaded regions in (a) indicate where the base state is asymptotically stable for $m=1000$.