2.1. Governing equations for the full 2-D model

The dimensionless governing equations for the fluid are the incompressible Navier–Stokes equations in the form

(2.6a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol u}&=0, \end{align}

(2.6b) \begin{align} \frac {\partial {\boldsymbol u}}{\partial t}+({\boldsymbol u}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol u}&=\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol \sigma }. \end{align}

The corresponding Newtonian stress tensor $\boldsymbol \sigma$ in (2.6b ) takes the form

(2.6c) \begin{align} {\boldsymbol \sigma }&=-p{\kern-1pt\boldsymbol{I}}+\textit{Re}^{-1}\left (\boldsymbol{\nabla}\kern-1pt{\boldsymbol u}+{\boldsymbol{\nabla}\kern-1pt{\boldsymbol{u}}}^{\text T}\right)\!, \end{align}

where $\boldsymbol I$ is the identity tensor and the superscript $^{\text T}$ denotes the transpose operator.

2.2. Boundary conditions

At a given time $t$ the flexible wall profile is the line of material points

(2.6d) \begin{equation} \{ (x_b(l,t),y_b(l,t)): 0\leqslant l \leqslant L \}. \end{equation}

Our numerical method (see § 2.4) can handle highly deformed beam configurations, including cases where the wall profile is a multi-valued function of $x$ (see the discussion in Luo & Pedley Reference Luo and Pedley1996). However, for all the examples considered below, it emerges that the $x$ component of the beam profile, $x_b$ , is a one-to-one function of $l$ for a given $t$ , so the dependency on $l$ can be eliminated and the corresponding profile of the wall can (in this case) be written in the form $y=h(x,t)=y_b(x_b^{-1}(x,t),t)$ .

For the viscous fluid, we assume no-slip conditions along the rigid walls and continuity of velocity between the fluid and elastic beam in the form

(2.6e) \begin{align} {\boldsymbol u}&={\boldsymbol 0}, & &\left (y=0;\,\,y=1,\,\, -L_u\leqslant x\leq 0, L\leqslant x\leqslant L+L_d\right ), \end{align}

(2.6f) \begin{align} {\boldsymbol u}&={\boldsymbol u}_s, & &\left (y=h(x,t),\,\,0\leqslant x \leqslant L\right ). \end{align}

For notational convenience, we denote $h(x,t)=1$ along $-L_u\leqslant x\leqslant 0$ and $L \leqslant x \leqslant{} L+L_d$ .

The governing equations for the massless elastic beam are obtained by considering conservation of linear and angular momentum (Cai & Luo Reference Cai and Luo2003; Wang et al. Reference Wang, Luo and Stewart2021a ). In turn, these can be expressed in the form of a balance of tangential and normal stress across the beam in the form

(2.6g) \begin{align} c_{\kappa }\frac {\kappa }{\lambda }\frac {\partial (\lambda \kappa )}{\partial l}+c_{\lambda }\frac {1}{\lambda }\frac {\partial \lambda }{\partial l}+\left ({\boldsymbol \sigma }\kern-1pt{\boldsymbol n}\right )\boldsymbol{\cdot }{\boldsymbol t}&=0, \end{align}

(2.6h) \begin{align} -{c_{\kappa }}\frac {1}{\lambda }\frac {\partial }{\partial l}\left (\frac {1}{\lambda }\frac {\partial (\lambda \kappa )}{\partial l}\right )+c_{\lambda }\kappa (\lambda -1)+\kappa T+\left ({\boldsymbol \sigma }\kern-1pt{\boldsymbol n}\right )\boldsymbol{\cdot }{\boldsymbol n}&=p_e, \end{align}

where we have the constraints

(2.6i) \begin{equation} \frac {\partial x_b}{\partial l}=\lambda \cos \theta ,\qquad \frac {\partial y_b}{\partial l}=\lambda \sin \theta ,\qquad \frac {\partial \theta }{\partial l}=\lambda \kappa ; \end{equation}

here $\lambda$ is the principal stretch of the elastic beam, $\kappa$ is the in-plane curvature and $\theta$ is the angle between the rigid wall and the tangent vector of the deformed elastic wall (see figure 1). The normal and tangential unit vectors of the deformed elastic beam are denoted as $\boldsymbol n$ and $\boldsymbol t$ , respectively.

In this paper we focus on a parameter regime where the elastic response of the beam is dominated by the longitudinal pre-tension, where the normal stress balance (2.6h ) can be approximated as

(2.7a) \begin{align} R \equiv C_T + C_n - p_e \approx 0, \end{align}

where we define the contributions due to the pre-tension and the normal fluid stress as

(2.7b) \begin{equation} C_T\equiv \kappa T, \qquad C_n\equiv \left ({\boldsymbol \sigma }\kern-1pt{\boldsymbol n}\right )\boldsymbol{\cdot }{\boldsymbol n}. \end{equation}

We quantify the residual $R$ and the size of these individual terms in figure 2.

Figure 2. Steady solutions of the fluid-beam system computed using the full 2-D simulations, plotting (a) the maximal and minimal steady wall deflection as a function of Reynolds number for $p_e=1$ ; (b) the spatial profile of the steady wall as a function of the material coordinate $l$ for $p_e=1$ , $Re=300$ ; (c) terms in the reduced normal stress balance (2.7a ) as a function of the streamwise coordinate $x$ for $p_e=1$ , $Re=300$ ; (d) the maximal and minimal steady wall deflection as a function of Reynolds number for $p_e=5$ ; (e) the spatial profile of the steady wall as a function of the material coordinate $l$ for $p_e=5$ , $Re=90$ ; ( f) terms in the reduced normal stress balance (2.7a ) as a function of the streamwise coordinate $x$ for $p_e=5$ , $Re=90$ . The upper and lower branch limit points for $p_e=1$ are denoted as filled squares in (a), while the approximate transition between the upper and lower branch for $p_e=5$ is indicated as an open square in (d). The point denoted with a circle in (a) corresponds to the profiles in (b,c), while the point denoted with a circle in (d) corresponds to the profiles in (e, f).

We note that in this study we have neglected the contribution of the wall inertia in the stress balances in the beam and so this prevents flutter-like instabilities found in previous studies (e.g. Luo & Pedley Reference Luo and Pedley1998) that onset due to two-way fluid–structure interaction. In addition, the two ends of the elastic beam are fixed and clamped to the rigid segments upstream and downstream, where we have

(2.8a) \begin{align} x_b(0,t)&=0, & y_b(0,t)&=1, & \theta (0,t)&=0, \end{align}

(2.8b) \begin{align} x_b(L,t)&=L, & y_b(L,t)&=1, & \theta (L,t)&=0. \end{align}

In this fluid-driven system the inlet flow has a prescribed parabolic profile $u({-L_u},y,t) = 6 y (1-y)$ , $v(-L_u,y,t)=0$ , with maximal velocity $3/2$ along the centreline. The corresponding fluid pressure along the channel inlet is denoted as $p_{\textit{in}}(y,t)=p(-L_u,y,t)$ ; in simulations we track temporal variations in the cross-sectionally averaged upstream pressure $\overline {p}_{\textit{in}}(t)$ (see details below).

In the computations we impose that the flow across the channel outlet ( $x=L+L_d$ ) is stress free, which for a sufficiently long downstream segment, is indistinguishable from a zero pressure condition (see discussion in Jensen & Heil Reference Jensen and Heil2003).

In the simulations reported below we track the dynamics of the unsteady system using several spatial and temporal measures. Firstly, we introduce the cross-sectionally averaged fluid pressure (denoted with an overbar), defined as

(2.9a) \begin{align} \overline {p}(x,t) =\frac {1}{h(x,t)}\int _0^{h(x,t)} p(x,y,t)\,\textrm {d}y & &(-L_u \leqslant x \leqslant L+L_d). \end{align}

In particular, across the channel inlet $\overline {p}_{\textit{in}}(t)=\overline {p}(-L_u,t)$ . Secondly, we consider the time trace of the wall pressure at the mid-point of the compliant segment in the reference configuration ( $l=(1/2)L$ ), which we denote as

(2.9b) \begin{equation} p_{\textit{mid}}(t) = p\left(x_b\left(\frac{1}{2}L,t\right), y_b\left(\frac{1}{2}L,t\right),t\right). \end{equation}

Thirdly, we consider the maximal and minimal elastic wall deflections denoted as

(2.9c) \begin{equation} y_{\textit{max}}(t)=\max _l(y_b(l,t))\geqslant 1, \qquad y_{\textit{min}}(t) = \min _l(y_b(l,t))\leqslant 1. \end{equation}

Additionally, the spatially averaged wall deflection over the compliant segment is denoted as

(2.9d) \begin{align} \overline {y}_b=\frac {1}{L}\int _0^L y_b(l,t)\,\textrm {d}l. \end{align}

Finally, we use the cross-sectionally averaged pressure distribution (2.9a ) to construct the spatially averaged pressure gradient along each compartment of the channel, such that

(2.9e) \begin{align} \overline {\boldsymbol{\nabla }\!p}_u(t) &= \frac {1}{L_u}\int _{-L_u}^0\ \frac {\partial\kern-1pt\overline{p}(x,t)}{\partial x} \,\textrm {d}x = \frac {1}{L_u}\left ({\overline {p}(0,t)-\overline {p}(-L_u,t)}\right )\!,\\[-8pt] \nonumber \end{align}

(2.9f) \begin{align} \overline {\boldsymbol{\nabla }\!p}_c(t) &= \frac {1}{L}\int _{0}^{L}\frac {\partial\kern-1pt\overline{p}(x,t)}{\partial x}\,\textrm {d}x = \frac {1}{L}\left ({\overline {p}(L,t)-\overline {p}(0,t)}\right )\!, \\[-8pt] \nonumber \end{align}

(2.9g) \begin{align} \overline {\boldsymbol{\nabla }\!p}_d(t) &= \frac {1}{L_d}\int _{L}^{L+L_d}\frac {\partial\kern-1pt\overline{p}(x,t)}{\partial x} \,\textrm {d}x = \frac {1}{L_d}\left ({\overline {p}(L_d,t)-\overline {p}(L,t)}\right )\!. \end{align}

For each of these temporal measures, we can then further compute the corresponding time-averaged value over one period of oscillation, which we denote with a superscript $^{(\textit{a}v\textit{g})}$ ; a function $f(t)$ averaged over a period $\tau$ takes the form

(2.10a) \begin{equation} f^{(\textit{a}v\textit{g})} = \frac {1}{\tau }\int _{t_0}^{t_0+\tau } f\kern-1pt(t)\,\textrm {d}t, \end{equation}

where $t_0$ is some reference time.

In our previous work (Wang et al. Reference Wang, Luo and Stewart2021a ) we also quantified the dynamics of self-excited oscillations using analysis of the corresponding energy budget of the system. We do not reproduce this derivation here but note that the rate of working of the upstream pressure can be expressed as

(2.10b) \begin{align} \mathcal{P}(t)=-\left [\int _0^1pu\,\textrm {d}y\right ]^{x=L+L_d}_{x=-L_u}. \end{align}

In the particular case where the cross-channel profile of the upstream driving pressure $p(-L_u,y,t)$ is independent of $y$ (which is a good approximation for the examples illustrated below), the rate of working of the upstream pressure $\mathcal{P}(t)\approx \overline {p}_{\textit{in}}(t)$ since the inlet flow rate is fixed and the outlet pressure is zero. Hence, the work done by the upstream pressure over a period of oscillation is given by

(2.10c) \begin{align} \mathcal{P}^{(\textit{a}v\textit{g})}\approx \overline {p}_{\textit{in}}^{(\textit{a}v\textit{g})}. \end{align}

2.3. Steady solutions

We denote the steady flow velocity, pressure and Newtonian stress tensor as $\boldsymbol{U}(x,y)$ , $P(x,y)$ and $\boldsymbol{\varSigma }(x,y)=-P(x,y){\boldsymbol I}+\textit{Re}^{-1} (\boldsymbol{\nabla }{\boldsymbol U}(x,y)+\boldsymbol{\nabla }{{\boldsymbol U}(x,y)}^{\text{T}} )$ respectively. The corresponding stretch and curvature of the steady beam are denoted as $\varLambda (l)$ and $K(l)$ , respectively, while the steady beam deflection is denoted as ${\boldsymbol X}_b(l)=(X_b(l), Y_b(l))$ and the angle between the steady beam and rigid wall is denoted as $\varTheta (l)$ . Therefore, the equations governing the steady system are

(2.11a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol U}&=0, \end{align}

(2.11b) \begin{align} \left ({\boldsymbol U}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ){\boldsymbol U}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\big(-P\kern-1pt{\boldsymbol I}+\textit{Re}^{-1}\left (\boldsymbol{\nabla }{\boldsymbol U}+\boldsymbol{\nabla }{\boldsymbol U}^{\text{T}}\right )\big), \end{align}

subject to boundary conditions

(2.11c) \begin{align} {\boldsymbol U}={\boldsymbol 0}& &\left (y=0;\,\,y=1,\,\, -L_u\leqslant x\leq 0, L\leqslant x\leqslant L+L_d; y=H(x),\,\,0\leqslant x\leqslant L\right )\!, \end{align}

and the governing equations for the steady beam are

(2.11d) \begin{align} c_{\kappa }\frac {K}{\varLambda }\frac {\partial\kern-1pt \left (\kern-1.5pt\varLambda\kern-1pt K\right )}{\partial l}+c_{\lambda }\frac {1}{\varLambda }\frac {\partial\kern-1pt\varLambda }{\partial l}+\left ({\boldsymbol \varSigma }\kern-1pt{\boldsymbol N}\right )\boldsymbol{\cdot }{\boldsymbol T}&=0, \end{align}

(2.11e) \begin{align} -c_{\kappa }\frac {1}{\varLambda }\frac {\partial }{\partial l}\left (\frac {1}{\varLambda }\frac {\partial\kern-1pt \left (\kern-1.5pt\varLambda\kern-1pt K\right )}{\partial l}\right )+c_{\lambda }K(\varLambda -1) + \textit{KT}+\left ({\boldsymbol{\varSigma }}{\boldsymbol N}\right )\boldsymbol{\cdot }{\boldsymbol N}&=p_e, \end{align}

where $\boldsymbol T$ and $\boldsymbol N$ are the corresponding tangential and normal unit vectors to the steady beam. In addition, the geometric constraints on the steady beam take the form

(2.11f) \begin{equation} \frac {\partial\kern-1pt X_b}{\partial l}=\lambda \cos \varTheta ,\qquad \frac {\partial\kern-1pt Y_b}{\partial l}=\varLambda \sin \varTheta ,\qquad \frac {\partial \varTheta }{\partial l}=\varLambda\kern-1pt K. \end{equation}

Similar to the unsteady system, we construct the steady wall profile in the form $y=H\kern-0.5pt(x)$ , and define $H\kern-0.5pt(x)=1$ along $-L_u \leqslant x \leqslant 0$ and $L \leqslant x \leqslant L+L_d$ for notational convenience. We form the steady cross-stream averaged pressure distribution

(2.12a) \begin{align} \overline {P}(x) = \frac {1}{H(x)}\int _0^{H(x)} P(x,y)\,\textrm {d}y & &(-L_u \leqslant x \leqslant L+L_d), \end{align}

alongside the steady inlet pressure $\overline {P}_{\textit{in}}=\overline {P}(-L_u)$ . Similarly, the spatially averaged steady pressure gradient along each compartment can be written as

(2.12b) \begin{align} \overline{\boldsymbol{\nabla}\!P}_u &= \frac {1}{L_u}({\overline {P}(0)-\overline {P}(-L_u)}), \end{align}

(2.12c) \begin{align} \overline{\boldsymbol{\nabla}\!P}_c &= \frac {1}{L}({\overline {P}(L)-\overline {P}(0)}), \end{align}

(2.12d) \begin{align} \overline{\boldsymbol{\nabla}\!P}_d &= \frac {1}{L_d}({\overline {P}(L_d)-\overline {P}(L)}). \end{align}

We also construct the maximal and minimal channel widths, defined as

(2.12e) \begin{equation} Y_{\textit{max}}=\max _l(Y_b(l))\geqslant 1, \qquad Y_{\textit{min}} = \min _l(Y_b(l))\leqslant 1, \end{equation}

respectively.

Finally, at a given time $t$ we compute the difference between the current wall profile and the steady profile, written as

(2.12f) \begin{equation} \Delta y_b(l,t) = y_b(l,t) - Y_b(l), \end{equation}

the difference between the steady and unsteady cross-stream averaged pressure distributions in the form

(2.12g) \begin{equation} \overline {\Delta p}(x,t) = \overline {p}(x,t)-\overline {P}(x), \end{equation}

as well as the excess (cross-sectionally averaged) streamwise pressure gradients along each compartment

(2.12h) \begin{align} \overline {{\delta p}}_{u}(t)&=\overline {\boldsymbol{\nabla }\!p}_u(t) -\overline{\boldsymbol{\nabla}\!P}_u, \end{align}

(2.12i) \begin{align} \overline {{\delta p}}_{c}(t) &= \overline {\boldsymbol{\nabla }\!p}_c(t) -\overline{\boldsymbol{\nabla}\!P}_c, \end{align}

(2.12j) \begin{align} \overline {{\delta p}}_{d}(t)&=\overline {\boldsymbol{\nabla }\!p}_d(t) -\overline{\boldsymbol{\nabla}\!P}_d. \end{align}

2.4. Numerical methods

Following Wang et al. (Reference Wang, Luo and Stewart2021a ) (see also Luo et al. Reference Luo, Cai, Li and Pedley2008) we use a finite element method with a six-node triangular mesh to solve the fluid-beam system (2.6). In particular, the flow is described using an Arbitrary Lagrangian–Eulerian method for the domain under the elastic wall with an adaptive mesh with rotational spines, whereas the domain under the rigid wall is described using an Eulerian method with a fixed mesh (Luo et al. Reference Luo, Cai, Li and Pedley2008). We employ a Galerkin approach to construct weighting functions for the fluid momentum equations (2.6a , 2.6b ) and the normal stress balance in the beam (2.6h ), but instead use a Petrov–Galerkin weighted residual method for the tangential stress balance across the beam (2.6g ) and the constraints (2.6i ). For unsteady problems, an implicit finite-difference scheme is used for time stepping, where a frontal approach and a Newton–Raphson iteration scheme are used at each time step to assemble and solve the discretised matrix equations.

To initialise the unsteady numerical solutions, we apply a small perturbation by using the corresponding steady state with a $1\,\%$ increment in $c_{\lambda }$ .

A mesh of 36 657 six-node triangular elements and 140 three-node elements is used for the numerical simulations. Convergence tests were reported in Wang et al. (Reference Wang, Luo and Stewart2021a ).

2.5. Eigenvalue problem

To further our understanding of the system, we also construct a linear stability eigensolver around a given steady state of the system. This eigensolver is constructed using the method described by Hao et al. (Reference Hao, Cai, Roper and Luo2016), where we add a small perturbation to the basic state in the form

(2.13) \begin{align} f(x,y,t)=F(x,y)+\varepsilon\kern-1pt\hat f(x,y,t),\quad g(l,t)=G(l)+\varepsilon\kern-1pt\hat g(l,t), \end{align}

where $\varepsilon \ll 1$ , $f\kern-0.5pt(x,y,t)$ and $g(l,t)$ represent the fluid and beam related variables, while $F\kern-0.5pt(x,y)$ , $G(l)$ and $\hat f(x,y,t)$ , $\hat g(l,t)$ are the corresponding steady state and first-order perturbation variables, respectively. To derive the corresponding eigenvalue problem for this perturbed system, we neglect all terms of $O(\epsilon ^2)$ and assume that these first-order perturbations take the wave-like form

(2.14) \begin{align} \hat f(x,y,t)=\tilde f(x,y)\textrm{e}^{-i\omega t}+\big(\tilde f(x,y)\textrm{e}^{-i\omega t}\big)^{*},\quad \hat g(l,t)=\tilde g(l)\textrm{e}^{-i\omega t}+\big(\tilde g(l)\textrm{e}^{-i\omega t}\big)^*, \end{align}

where $\tilde f(x,y)$ and $\tilde g(l)$ are the complex amplitude functions and $\omega =\omega _r+i\omega _i$ is the complex frequency; here $\omega _r$ is the oscillation frequency of the perturbation and $\omega _i$ is the corresponding growth rate. The system is unstable when $\omega _i\gt 0$ , asymptotically stable when $\omega _i\lt 0$ and is neutrally stable for $\omega _i=0$ . The superscript $^*$ denotes a complex conjugate.

Substituting (2.13) and (2.14) into the governing equations (2.6), we obtain the eigenvalue problem for the coupled system in the form

(2.15a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol u}&=0, \end{align}

(2.15b) \begin{align} -i\omega \tilde {\boldsymbol u}+\left ({\boldsymbol U}\boldsymbol{\cdot }\boldsymbol{\nabla }\right )\tilde {\boldsymbol u}+\left (\tilde {\boldsymbol u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ){\boldsymbol U}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\left[-\tilde p\kern-0.5pt{\boldsymbol I}+\textit{Re}^{-1}\left (\boldsymbol{\nabla }\tilde {\boldsymbol u}+\boldsymbol{\nabla }\tilde {\boldsymbol u}^{\text T}\right )\right ]\!, \end{align}

subject to boundary conditions (projected to the steady state)

(2.15c) \begin{align} &\tilde {\boldsymbol u}={\boldsymbol 0} & &\left (y=0;y=1, -L_u\leqslant x\leq 0, L\leqslant x\leqslant L+L_d\right )\!, \end{align}

(2.15d) \begin{align} &\tilde {\boldsymbol u}+\left (\tilde {x}_b\frac {\partial {\boldsymbol U}}{\partial x}+\tilde {y}_b\frac {\partial {\boldsymbol U}}{\partial y}\right )=\tilde {\boldsymbol u}_s(x) & &\left (y=H, 0\leqslant x\leqslant L\right )\!, \end{align}

where $\tilde {\boldsymbol u}_s(x)=-i\omega \tilde {\boldsymbol x}_b ({\boldsymbol X}^{-1}_b(x) )$ is the linearised spatial wall velocity. The linearised governing equations for the elastic beam take the form

(2.15e) \begin{align} &c_{\kappa }\left [\frac {\tilde \kappa }{\varLambda }\frac {\textrm {d}(\varLambda K)}{\textrm {d} l}+\frac {K}{\varLambda }\frac {\textrm {d}(\tilde \lambda K+\varLambda \tilde \kappa )}{\textrm {d} l}-\frac {\tilde \lambda K}{{\varLambda }^2}\frac {\textrm {d}(\varLambda K)}{\textrm {d}l}\right ]+c_{\lambda }\left (\frac {1}{\varLambda }\frac {\textrm {d}\tilde \lambda }{\textrm {d} l}-\frac {\tilde \lambda }{{\varLambda }^2}\frac {\textrm {d}\varLambda }{\textrm {d}l}\right )\notag \\ &\qquad +\left [\left (\tilde {\boldsymbol{\sigma }}+\tilde {x}_b\frac {\textrm {d}{\boldsymbol{\varSigma }}}{\textrm {d} x}+\tilde {y}_b\frac {\textrm {d}{\boldsymbol \varSigma }}{\textrm {d}y}\right ){\boldsymbol N}\right ]\boldsymbol{\cdot }{\boldsymbol T}=0, \end{align}

(2.15f) \begin{align} &-c_{\kappa }\left [\frac {1}{\varLambda }\frac {\textrm {d}}{\textrm {d}l}\left (\frac {1}{\varLambda }\frac {\textrm {d}(\varLambda \tilde \kappa +\tilde \lambda K)}{\textrm {d}l}-\frac {\tilde \lambda }{{\varLambda }^2}\frac {\textrm {d}(\varLambda K)}{\textrm {d}l}\right )-\frac {\tilde \lambda }{{\varLambda }^2}\frac {\textrm {d}}{\textrm {d}l}(\frac {1}{\varLambda }\frac {\textrm {d}\left (\varLambda K\right )}{\textrm {d}l})\right ]\notag \\ &\qquad +c_{\lambda }\left [ K\tilde \lambda +\tilde \kappa \left (\varLambda -1\right )\right ]+\tilde \kappa T+\left [\left (\tilde {\boldsymbol{\sigma }}+\tilde {x}_b\frac {\textrm {d}{\boldsymbol{\varSigma }}}{\textrm {d} x}+\tilde {y}_b\frac {\textrm {d}{\boldsymbol \varSigma }}{\textrm {d}y}\right ){\boldsymbol N}\right ]\boldsymbol{\cdot }{\boldsymbol N}=0, \end{align}

subject to the constraints

(2.15g) \begin{align} &\frac {\textrm {d}\tilde {x}_b}{\textrm {d}l}=-\varLambda \tilde \theta \sin (\varTheta )+\tilde \lambda \cos (\varTheta ),\quad \frac {\textrm {d}\tilde {y}_b}{\textrm {d}l}=\varLambda \tilde \theta \cos (\varTheta )+\tilde \lambda \sin (\varTheta ),\quad \frac {\textrm {d}\tilde \theta }{\textrm {d}l}=\varLambda \tilde \kappa +\tilde \lambda K. \end{align}

We note that these perturbation eigenfunctions are all projected back to the corresponding steady state, which leads to a number of additional terms in the linearised governing equations involving spatial gradients of the steady flow vector $\boldsymbol{U}$ and the steady stress tensor $\boldsymbol{\varSigma }$ ; see Wang (Reference Wang2019) for details. We note that although these extra terms are crucial to the mechanism of instability in many classical hydrodynamic stability problems (e.g. Kelvin–Helmholtz and Rayleigh–Taylor instabilities; see Drazin Reference Drazin2002), in this case their influence is relatively small.

Following Hao et al. (Reference Hao, Cai, Roper and Luo2016), we solve the eigenvalue problem (2.15) numerically using an Arnoldi-Frontal method. We characterise the output of the eigensolver through the growth rate $\omega _i$ and frequency $\omega _r$ of the normal modes, alongside the corresponding spatial profiles of the eigenfunction (which has both real and imaginary parts).

2.6. Reduced 1-D model

We also consider a reduced 1-D model of the same collapsible channel system, first derived by Stewart et al. (Reference Stewart, Waters and Jensen2009) and subsequently used by a number of authors (e.g. Stewart et al. Reference Stewart, Heil, Waters and Jensen2010; Xu et al. Reference Xu, Billingham and Jensen2013, Reference Xu, Billingham and Jensen2014; Xu & Jensen Reference Xu and Jensen2015; Stewart Reference Stewart2017). The derivation of this model assumes long-wavelength deformations, a von Karman–Pohlhausen approximation for the flow profile (i.e. everywhere parabolic), encodes the reduced form of the normal stress balance (2.6h ) and neglects axial stretching of the wall (so motion of the wall is always normal to the rigid wall). In this model the system is described in terms of the width of the collapsible segment $h(x,t)$ and the axial flux along the channel $q(x,t)$ , where the governing equations can be written in the form

(2.16a) \begin{align} h_t + q_x =0,\quad q_t + \frac {6}{5}\left (\frac {q^2}{h} \right )_x = T hh_{xxx} - \frac {12}{\textit{Re}} \frac {q}{h^2}, \end{align}

across $0\leqslant x \leqslant L$ , subject to

(2.16b) \begin{equation} q(0,t) = 1, \qquad Th_{xx}(L,t) = - L_d \left [q_t + \frac {12}{\textit{Re}}{q}\right ]\!,\qquad h(0,t)=h(L,t)=1. \end{equation}

This system is solved numerically using finite-difference methods constructed by Stewart et al. (Reference Stewart, Waters and Jensen2009), using Newton iteration to construct the steady states and a semi-implicit time-stepping scheme to compute unsteady profiles. Furthermore, we use an analogue of our linear stability eigensolver (constructed by Stewart Reference Stewart2017) to test the stability of steady states, constructing the growth rate and frequency of perturbations, as well as their corresponding complex eigenfunctions.

2.7. Parameter choices

Throughout this study we fix the lengths of the three segments of the channel ( $L_u$ , $L$ and $L_d$ ) and the initial thickness of the elastic beam $\mathcal{H}$ by setting

(2.17) \begin{equation} L_u=L=5,\quad L_d=30, \quad \mathcal{H}=0.01. \end{equation}

This choice of downstream segment is sufficiently long so that outflow is unidirectional.

We focus attention on a parameter regime where the pre-tension in the wall strongly dominates the other elastic effects, meaning that the system is analogous to earlier membrane models (e.g. Luo & Pedley Reference Luo and Pedley1996, Reference Luo and Pedley2000). In this case we choose $c_{\lambda }=10$ , $c_{\kappa }=(\mathcal{H}^2/12)c_\lambda$ and $T=3$ ; we confirm in figure 2 that the dominant elastic effect is indeed the restoring effect of the pre-tension for these parameters.

To examine the steady and unsteady behaviour of the system, we choose the external pressure in the range $0 \lt p_e \leqslant 10$ and focus on laminar (non-turbulent) flows in the range $0 \lt \textit{Re} \lt 800$ .

In what follows we summarise the predictions of our models, including a brief overview of both the steady (see § 3) and unsteady behaviour (see § 4) and a more detailed study of particular families of oscillation (see §§ 5, 6, Appendix).

4.1. Neutral stability curves

Figure 3 presents an overview of the stability of the collapsible channel system, plotting the neutral stability curves in the parameter space spanned by external pressure $p_e$ and Reynolds number $\textit{Re}$ , holding the wall pre-tension fixed at $T=3$ ; the full space is shown in figure 3(a), while panels (b,c) provide enlargements of particular regions (magenta dashed boxes in figure 3 a). Tracing in terms of the external pressure, the corresponding maximal ( $Y_{\textit{max}}$ ) and minimal ( $Y_{\textit{min}}$ ) channel widths of the steady profile at neutral stability are plotted in figure 3(d), while the neutrally stable oscillation frequencies are plotted in figure 3(e).

For large external pressures, the system exhibits a neutral stability curve where the critical Reynolds number decreases with increasing $p_e$ (see figure 3 a), where the corresponding steady channel is highly collapsed (consistent with the lower branch of steady solutions, figure 3 d) and the frequency of the neutrally stable oscillation is $O(1)$ (figure 3 e). The corresponding nonlinear growth of these lower branch oscillations is discussed in § 5. The onset of this lower branch oscillatory instability is also well predicted by the corresponding 1-D model (2.16), where the neutral stability curve, static profiles and oscillation frequencies all agree well with the full computations (figure 3 d,e). Furthermore, we also see strong agreement with the asymptotic solution of this 1-D model in the limit of large $p_e$ (Stewart Reference Stewart2017), where the critical Reynolds number for the onset of oscillatory instability takes the approximate form

(4.1) \begin{equation} Re_c \approx \frac {12}{p_e}L_d; \end{equation}

this curve is shown as the black dot-dashed line in figure 3(a), which overlaps closely with the other approaches. Hence, the 1-D model is a strong predictor of the oscillatory behaviour along this lower branch of steady solutions.

As the external pressure decreases, the computed neutral stability curve becomes increasingly complicated. The lower branch oscillatory structure is eventually suppressed at a limit point ( $p_e \approx$ 1.80), beyond which the system transitions to upper branch steady states. In this case the limit point is one of a pair (the other being at $p_e \approx$ 2.55) that form a sigmoid in the neutral stability curve across a short range of external pressures ( $1.80\lesssim p_e \lesssim 2.55$ ), where the system exhibits three possible critical Reynolds numbers (see close-up of the parameter space in figure 3 c). Across these three neutrally stable solutions, the corresponding steady profile becomes increasingly inflated as the Reynolds number decreases (figure 3 d), while the corresponding oscillation frequency gradually increases (figure 3 e). We explore the nonlinear development of several oscillatory modes in this regime in § 4.2. We note that these limit points are not associated with those of the underlying steady system (those limit points only arise for much larger Reynolds numbers, shown as the shaded narrow red region in figure 3 a), but our approximate demarcation line between upper and lower branches (red dashed line) does pass almost perfectly through the second limit point at $p_e\approx 2.55$ , suggesting that this point is indeed associated with a change from a lower to an upper branch steady configuration. Note also that the 1-D model (2.16) does not predict this three branch behaviour, and instead the corresponding neutral stability curve persists in the lower branch configuration until eventually terminating when the system reaches a limit point associated with the end of the lower steady branch (this point is indicated with a filled blue square in figures 3 a and 3 b).

Returning to the neutral stability curve from the full computations, we observe that as the external pressure decreases below the lower limit point ( $p_e \approx 2.55$ ), the critical Reynolds number again approximately increases with the inverse of $p_e$ , although the corresponding steady profile is almost entirely inflated and so the dynamics of the self-excited oscillations growing from this upper branch are quite different to the lower branch. We consider the nonlinear development of oscillatory normal modes growing from this upper branch in § 6.

However, across the range of external pressures with three neutrally stable solutions we detect a second, entirely separate, neutral stability curve associated with an instability of much higher frequency growing from the lower branch of steady solutions (marked as the green line in figure 3). Across a narrow range of external pressures ( $1.5 \lesssim p_e \lesssim 3.3$ ) this high-frequency normal mode is the primary global instability, i.e. more unstable than the low-frequency mode. A narrow range of the external pressures where the collapsible channel system is dominated by high-frequency normal modes was also evident in the predictions from the 1-D model (Stewart Reference Stewart2017). This high-frequency family of modes is described in more detail in § 5.2.

For this low value of the membrane tension, it emerges that the upper branch instability is more complicated than previously reported by Wang et al. (Reference Wang, Luo and Stewart2021a ), where the instability was eventually suppressed as the external pressure decreased below a threshold (Wang et al. Reference Wang, Luo and Stewart2021b ). In the present case the system exhibits a second form of upper branch instability for even lower values of the external pressure, where the dynamics of the resulting self-excited oscillations are substantially different. The corresponding neutral stability curve for this instability is evident in the parameter space in figure 3, manifesting as a large tongue (coloured in purple) sitting below the neutral stability curve discussed above, connecting smoothly to it at $p_e\approx 0.25, Re\approx 805.82$ . Furthermore, this new two-branch neutral stability curve is confined to low external pressures by a limit point at $p_e \approx 4.02$ , beyond which this new form of upper branch instability is absent. We consider an oscillatory normal mode in this unstable range in the Appendix.

4.2. Traces through the parameter space

The parameter space shown in figure 3 demonstrates that the system is unstable to at least three distinct forms of self-excited oscillation (termed lower branch, upper branch and fully upper branch oscillations).

In order to show how these different families of instabilities fit together with their corresponding static profiles, figure 4 presents several slices through the parameter space for a fixed external pressure. In particular, we construct bifurcation diagrams as a function of Reynolds number by plotting the time-averaged mid-point pressure on the compliant wall, $p_{\textit{mid}}^{(\textit{a}v\textit{g})}$ , for three external pressures, namely $p_e=5$ (figure 4 a), $p_e=2.2$ (figure 4 b) and $p_e=1$ (figure 4 c) alongside the corresponding values computed from the steady simulations (where stable and unstable configurations are plotted as solid and dashed lines, respectively).

Figure 4. Overview of the parameter space calculated using full computations of the 2-D model, plotting both the steady and time-averaged mid-point pressures as a function of the Reynolds number for (a) $p_e=5$ , (b) $p_e=2.2$ and (c) $p_e=1$ . In each panel the time-averaged mid-point pressures for oscillations growing from the upper and lower steady branches are denoted as filled black and magenta circles, respectively, while the stable and unstable steady states are denoted with solid and dashed lines, respectively. The approximate demarcation between the upper and lower steady states is marked with an open red square in panel (b). The corresponding steady and time-averaged mid-point pressures calculated from the 1-D model (2.16) are denoted in blue with the same line styles.

For relatively large external pressures ( $p_e=5$ , figure 4 a), no upper branch or fully upper branch instability is present and the system becomes unstable to low-frequency lower branch oscillations beyond a critical Reynolds number ( $Re \approx 87.17$ , labelled N6). The corresponding transition point for the 1-D model is slightly lower ( $Re \approx 86.25$ ), but the neutrally stable eigenfunctions show excellent qualitative agreement (see figure 5). These oscillations grow in amplitude as the Reynolds number increases (figure 4 a) and eventually the system exhibits highly nonlinear ‘slamming’ oscillations for $Re \gtrsim 95$ for the full simulations (marked by pink arrow in figure 4 a). Fully nonlinear simulations of the 1-D model (2.16) exhibit qualitatively similar dynamics, although the growth in amplitude of the saturated limit cycle is much more pronounced as the Reynolds number increases, and so the ‘slamming’ oscillations onset much closer to the neutral stability curve ( $Re \gtrsim 87.5$ in this case).

Figure 5. Neutrally stable low-frequency self-excited oscillations originating from the lower branch of steady solutions for $p_e=5$ , considering the prediction of the 2-D model at point N6 ( $Re=87.17$ ) alongside the prediction from the 1-D model ( $Re=86.25$ ): (a) steady flow-field and pressure contours from the 2-D model, (b) eigenfunction wall profile plotted at 10 snapshots over a period of oscillation from the 2-D model, (c) eigenfunction wall profile plotted for 10 snapshots over a period of oscillation from the 1-D model. The eigenfunctions in (b,c) are normalised on the maximum of the modulus of the (complex-valued) wall profile. The blue and cyan filled circles in (a) are the corresponding cross-stream averaged steady pressure distributions for the one- and 2-D models, respectively. The dashed lines in (a–c) indicate the spatial location of the minimal channel width, while the dot-dashed lines in (a–c) indicate the point where the wave profile is approximately pinned; in each case the magenta lines are from the 2-D model, while the blue lines are from the 1-D model. The insets in (b,c) show the steady wall profile from 2-D and 1-D models, respectively.

For an external pressure chosen within the interval that exhibits three pieces of the same neutral stability curve ( $p_e=2.2$ , figure 4 b), the system now exhibits two families of oscillatory behaviour. For low Reynolds numbers, the system exhibits an unstable interval between $Re \approx 172.05$ (labelled N3) and $Re \approx 180.97$ (labelled N4). The corresponding steady profile of the channel is completely inflated for lower Reynolds numbers across this range, but becomes partially collapsed as the Reynolds number increases. These oscillations are similar in character to the upper branch oscillations described in our previous work (Wang et al. Reference Wang, Luo and Stewart2021a ), and while their eventual termination as the Reynolds number increases is not directly associated with the limit point of the steady solutions, it is driven by the collapse of the underlying steady state (similar to examples shown in Wang et al. Reference Wang, Luo and Stewart2021b ). Indeed this normal mode growing from the upper branch of steady solutions stabilises close to our approximate demarcation between upper and lower steady states (open red square in figure 4 b). The mechanism of these upper branch instabilities is explored in § 6. As the Reynolds number increases, the system subsequently destabilises again for $Re \approx 184.63$ (labelled N5); this transition is analogous to the onset of lower branch oscillations described for $p_e=5$ , but corresponds to a different normal mode of the system with a much larger oscillation frequency (figure 3 e); the temporal dynamics of these high-frequency lower branch oscillations are explored in § 5.2.

For a relatively low external pressure ( $p_e=1$ , figure 4 c), the system again exhibits two intervals of oscillatory behaviour. In this case the oscillatory instability first onsets for $Re \approx 250.01$ (labelled N1), where the corresponding steady state is very highly inflated (figure 3 d). As the Reynolds number increases, this oscillatory mode subsequently stabilises at $Re\approx 273.02$ , where the base state is still very highly inflated. Hence, the dynamics of these fully upper oscillations is different to the upper branch oscillations described previously and are explored in more detail in the Appendix. As the Reynolds number increases further, the system eventually destabilises again for $Re \approx 343.67$ (labelled N2), in a manner almost identical to the upper branch oscillations. For this external pressure, the system exhibits a range of Reynolds numbers with three coexisting steady states. For these parameters, the corresponding lower branch of steady solutions is always unstable to oscillations, and the associated limit cycle merges with that of the upper branch oscillations across the interval with multiple steady states; a similar merging of the upper and lower branch limit cycles was observed by Wang et al. (Reference Wang, Luo and Stewart2021b ).

We now proceed to explore the nonlinear dynamics of the lower branch oscillations (§ 5), the upper branch oscillations (§ 6) and the fully upper branch oscillations (Appendix) by considering 12 particular cases (labelled on figure 4), elucidating the underlying mechanism of instability by examining both the neutrally stable eigenfunctions (operating points N1–N6) and the unsteady dynamics nearby (operating points P1–P6). Movies of the unsteady oscillations at points P1–P6 are provided as supplementary material available at https://doi.org/10.1017/jfm.2025.10860.

5.1. Low-frequency normal modes

These lower branch oscillations are analogous to the low-frequency lower branch oscillations described by Stewart (Reference Stewart2017) using the (empirical) 1-D model (2.16). In order to validate these predictions against full 2-D computations, figure 5 explores the dynamics of neutrally stable low-frequency lower branch oscillations at operating point N6. Firstly, we plot streamlines and pressure contours of the steady flow profile and the cross-sectionally averaged pressure distribution (figure 5 a) alongside the corresponding steady wall profile from the 1-D model (shown in blue, plotted at the critical Reynolds number $Re \approx 86.25$ for the same external pressure). These steady wall profiles compare well between the two approaches, both exhibiting an outward bulge at the upstream end of the compliant compartment, and substantial collapse toward the downstream end in a so-called mode-2 configuration (two extrema in the wall profile), while the locations of the minimal channel thickness are also very similar (figure 5 a). The corresponding steady pressure profiles also show good agreement, both exhibiting a local minimum toward the downstream end of the compliant segment. Snapshots of the eigenfunction wall profile at 10 equally spaced points over a period of neutrally stable oscillation are shown for both the full simulations (figure 5 b) and the 1-D model (figure 5 c). Over a period of oscillation these neutrally stable eigenfunction profiles agree very well between the two approaches, exhibiting an almost identical amplitude (when normalised according to the same protocol), where each eigenfunction wall profile has two extrema (mode-2) with a larger amplitude toward the downstream end of the collapsible segment. Interestingly, in both approaches all the oscillatory wall profiles exhibit zero amplitude at a consistent spatial location just downstream of the centre of the compliant compartment (i.e. both real and imaginary parts of the eigenfunction wall profile are close to zero simultaneously, indicated by the dot-dashed lines in figures 5 b and 5 c); this zero-crossing point is distinct from the global minimum of the steady profile (which is considerably closer to the downstream end of the compliant channel, see dashed line in figures 5 a and 5 c), and so the steady state is not dictating the structure of the eigenfunction (in contrast to the high-frequency instabilities discussed in § 5.2). Although we note small variations in the location of this zero-crossing for 2-D simulations (figure 5 b), for the 1-D model, this point is almost perfectly constant (figure 5 c), suggesting that the perturbation wall profile is approximately a standing wave. In summary, this figure suggests that the neutrally stable perturbation wall profile is approximately a standing wave, pinned both at the ends of the compliant segment and also at an interior point near the centre.

In order to explore the dynamics of the fully developed limit cycles originating from this oscillatory normal mode, figure 6 illustrates temporal and spatial profiles of the flow over a period of self-excited oscillation for a Reynolds number just above neutral (operating point P6). A time trace of the mid-point wall pressure $p_{\textit{mid}}$ is shown in figure 6(a), illustrating three periods of oscillation, while figure 6(b–k) shows streamlines and pressure contours at the 10 points marked on this trace. For each snapshot, we also plot the corresponding spatial profile of the difference between the (cross-stream averaged) instantaneous pressure and the steady pressure (i.e. $\overline {\Delta p}(x,t)$ , figure 6 b–k). The flow profiles are almost perfectly laminar throughout the oscillation and exhibit little cross-stream pressure variation (the colours are almost independent of $y$ ), in agreement with the assumptions of the 1-D model (figure 6 b–k). In all cases the unsteady pressure gradient in the upstream rigid segment is essentially identical to the corresponding steady flow (the excess pressure profile is spatially uniform), while in the downstream rigid segment the unsteady pressure distribution varies almost linearly compared with the steady pressure distribution (the excess pressure profile is linear; see cyan triangles in figure 6 b–k). However, across the compliant segment variations in the pressure are more significant, in some cases almost half of the magnitude of the corresponding steady driving pressure ( $\overline {P}_{\!\textit{in}}$ ). We consider these quantities further in figure 7.

Figure 6. Fully developed low-frequency self-excited oscillations growing from the lower branch of steady solutions predicted by the 2-D model at operating point P6 ( $p_e=5$ , $Re=90$ ): (a) time trace of the wall mid-point pressure $p_{\textit{mid}}(t)$ over three periods of oscillation, (b–k) streamlines and pressure contours of the instantaneous flow field at 10 equally spaced time instances over one period (corresponding times marked in panel (a)). The cyan filled triangles in (b–k) illustrate the excess (cross-stream averaged) pressure compared with the steady pressure ( $\overline {\Delta}p$ ). The insets in (b–k) show the difference between the instantaneous wall profile compared with the steady flow ( $\Delta y_b$ ).

Figure 7. Fully developed low-frequency self-excited oscillations growing from the lower branch of steady solutions at operating point P6 ( $p_e=5$ , $Re=90$ ), including (a) the time trace of the cross-sectionally averaged inlet pressure ( $\overline {p}_{\textit{in}}$ , black curve) and the excess pressure gradient compared with the steady ( $\delta {p}_{\boldsymbol{\cdot }}$ ) in the upstream (blue curve), compliant (red curve) and downstream (green curve) segments; (b) a close-up of the time trace of the cross-sectionally averaged inlet pressure; (c) a close-up of the time trace of the excess pressure gradients in each segment; (d) the time trace of the minimal channel width (left axis) and the spatially averaged channel width (right axis). Note the corresponding time-averaged and steady values are denoted as dashed and dotted lines in panels (a–d).

Over an oscillation the wall retains the mode-2 wall profile with almost no longitudinal motion of the material points, while the location of maximal constriction also remains approximately fixed ( $x \approx 3.9$ ). As expected, at instances where the channel is highly collapsed the fluid pressure upstream of the point of maximal constriction becomes raised to become almost equal to the driving pressure (see colour contours in figure 6 e–h). Furthermore, as the channel collapses the pressure close to the point of maximal constriction falls considerably compared with the steady value (figure 6 d–g), but still remains well in excess of the outlet pressure and so there is no flow reversal in the downstream rigid segment. We also note that, for this relatively low Reynolds number, there is also no substantial vorticity wave downstream of the compliant segment, a prominent feature of earlier models in flux-driven oscillations in collapsible channel flows (e.g. Luo & Pedley Reference Luo and Pedley1996; Luo et al. Reference Luo, Cai, Li and Pedley2008; Huang et al. Reference Huang, Tian, Young and Lai2021). Hence, the resulting oscillation predicted by the full 2-D computations exhibits many of the features encoded into the 1-D model (2.16), and so it is no surprise that this 1-D model predicts quantitatively similar oscillations. The dynamics of the oscillation can be seen more clearly in the movie for operating point P6 provided as supplementary material.

Stewart (Reference Stewart2017) used an approximate 1-D model to describe the mechanism of these low-frequency lower branch oscillations, showing that on average over a period the compliant wall is less collapsed than the corresponding steady wall profile, which releases the energy needed to sustain the oscillation. In order to validate these predictions using the full 2-D system, figure 7 compares time traces of several pertinent quantities to their corresponding steady values, including the (cross-sectionally averaged) upstream driving pressure (figure 7 a,b), the (spatially averaged) streamwise pressure gradient along each segment of the channel (figure 7 a,c) and both the minimal width of the channel and the average wall position (figure 7 d). Over a period of oscillation the mean upstream pressure is reduced compared with the corresponding steady value, meaning that (on average) the oscillatory flow is working less hard than the steady flow, making it energetically favourable (figure 7 a,b). Furthermore, over a period of oscillation there is no significant change in the pressure gradient along the upstream and downstream rigid segments of the channel compared with the steady flow (figure 7 a,c). However, the corresponding time-averaged pressure gradient along the compliant segment is significantly increased (less negative, figure 7 a,c), corroborating our observation that overall the flow is working less hard in the oscillatory state. This reduction is also evident in the spatially averaged beam deflection, where, in agreement with Stewart (Reference Stewart2017), on average the beam profile is less collapsed over a period compared with the steady state and so stored elastic energy is released to sustain the oscillation (figure 7 d). Hence, this figure demonstrates that our fully nonlinear 2-D simulations agree well with the earlier predictions of the 1-D model, showing that on average the fully developed oscillation reduces the mean constriction of the channel and increases the pressure gradient along the compliant segment. In the online supplementary material we also demonstrate that the cross-channel profiles of the velocity components and the pressure are consistent with the assumptions of the 1-D model. In summary, figure 7 shows that the mechanism of these low-frequency self-excited oscillations growing from the lower branch of steady solutions is similar to that predicted by the 1-D model of Stewart (Reference Stewart2017), where the fully developed oscillation causes the collapsed wall to open slightly (on average), releasing stored elastic energy to sustain the instability.

5.2. High-frequency normal modes

In addition to these low-frequency normal modes, Stewart (Reference Stewart2017) also reported a family of high-frequency normal modes that grow from the lower branch of steady solutions. These normal modes are typically more stable than the low-frequency oscillations, but Stewart (Reference Stewart2017) elucidated a small region of the parameter space where one of these high-frequency modes is the primary global instability of the system. These oscillations were described using an asymptotic approximation in the limit of large external pressure, where it was shown that the steady channel becomes highly constricted, pinning the perturbation wall profile at the location of the global minimum of the steady profile, and the corresponding perturbation flow can thus be effectively divided into two pieces, each of which is characterised by the number of extrema in the perturbation wall profile. In order to compare these predictions to those of the full 2-D model, in figure 8 we consider the dynamics at operating point N5, a neutral stability point within the small range of parameter space where the high-frequency response is the primary global instability (see figure 3). As expected, the steady state is collapsed (figure 8 a), although the maximal constriction is relatively modest ( $Y_{\textit{min}}\approx 0.67$ at $X_{\textit{min}} \approx 3.54$ ). However, even in this case the corresponding perturbation wall profiles (figure 8 b) are standing waves pinned at the point of maximal constriction, with two (one) extrema (extremum) upstream (downstream) of the point of maximal constriction, in good qualitative agreement with one of the high-frequency oscillatory modes found using the 1-D model (Stewart Reference Stewart2017). Hence, in this case the perturbation wall profile exhibits two zero crossings with almost coincident real and imaginary parts (red dashed line in figure 8 b). In order to further elucidate the oscillatory dynamics, figure 8(c) traces the spatial location of the local minimum of the perturbation wall profile, demonstrating how it gradually transverses upstream following a path that takes two periods of oscillation. In summary, this figure demonstrates that the family of high-frequency oscillatory modes predicted by the 2-D simulations exhibit strong qualitative agreement with the earlier predictions of the 1-D model, where the wall profile is approximately a standing wave pinned at the point of maximal constriction of the steady profile.

Figure 8. Neutrally stable high-frequency self-excited oscillations originating from the lower branch of steady solutions for $p_e=2.2$ from the 2-D model at point N5 ( $Re=184.63$ ): (a) steady flow-field and pressure contours, (b) eigenfunction wall profile plotted at 10 snapshots over a period of oscillation from the 2-D model, (c) spatial location and value of the local minimum of the perturbation wall profile over a period of oscillation. The eigenfunctions in (b,c) are normalised on the maximum of the modulus of the (complex-valued) wall profile. The cyan filled circles in (a) are the corresponding cross-stream averaged steady pressure distributions. The red dashed lines in (a,b) indicate the spatial location of the minimal channel width, while the magenta dot-dashed line in (b) indicates the point where the wave profile is approximately pinned.

Figure 9. Fully developed high-frequency self-excited oscillations growing from the lower branch of steady solutions at operating point P5 ( $p_e=2.2$ , $Re=187$ ), including (a) the time trace of the cross-sectionally averaged inlet pressure ( $\overline {p}_{\textit{in}}$ , black curve) and the excess pressure gradient compared with the steady ( $\delta {p}_{\boldsymbol{\cdot }}$ ) in the upstream (blue curve), compliant (red curve) and downstream (green curve) segments; (b) a close-up of the time trace of the cross-sectionally averaged inlet pressure; (c) a close-up of the time trace of the excess pressure gradients in each segment; (d) the time trace of the minimal channel width (left axis) and the spatially averaged channel width (right axis). Note the corresponding time-averaged and steady values are denoted as dashed and dotted lines in panels (a–d).

For brevity we do not include snapshots of the streamlines and pressure contours, but the unsteady dynamics at operating point P5 are shown in the online supplementary material and in the accompanying movie. Interestingly, the excess cross-stream averaged pressure ( $\overline {\Delta p}(x,t)$ ) now typically exhibits three extrema over the compliant segment rather than two for the low-frequency oscillations (see figure 5).

In order to elucidate the mechanism of these high-frequency instabilities, in figure 9 we trace the (cross-sectionally averaged) upstream driving pressure (figure 9 a,b), the corresponding difference in the (spatially averaged) pressure gradient between unsteady and steady flows in each compartment (figure 9 a,c) as well as both the spatially averaged and minimal beam deflections (figure 9 d). In this case the effect is much more modest and harder to extract: the time-averaged upstream driving pressure is increased compared with the steady value (figure 9 b), while the other trends are similar to the lower branch oscillations in that the spatially averaged pressure is again increased compared with the steady profile (red dashed line in figure 9 a,c) and the time-averaged beam deflection is decreased compared with the steady value (figure 9 d). Hence, this figure demonstrates that the mechanism of the high-frequency lower branch instability is much more difficult to discern, exhibiting some overlap with the low-frequency lower branch oscillations but the observations are not fully consistent. However, given the relatively small differences between the steady and time-averaged unsteady profiles, more investigation is required to fully quantify the mechanism. In online supplementary material we again show that the cross-channel profiles of the velocity components and pressure are consistent with the assumptions of the 1-D model.