2.1 Static solutions

Static solutions to (2.8) must have unit flux to satisfy the inlet flow condition, while the corresponding membrane profile $h=h^{(s)}(x)$ must satisfy the third-order nonlinear ordinary differential equation (ODE)

subject to three boundary conditions

The superscript $(s)$ has been used throughout to denote the static solution and $p^{(s)}(x)$ denotes the corresponding static pressure profile. Following Stewart et al. (Reference Stewart, Waters and Jensen2009), these equations can be discretised using finite differences and solved numerically using Newton’s method. Static solutions to (2.9) are described in § 3. In a similar manner to Xu et al. (Reference Xu, Billingham and Jensen2013), it is also instructive to consider inviscid static solutions, discussed in appendix A. Throughout this paper we consider relatively large external pressures where the channel is collapsed for a region along its length, and we typically characterise the static solutions using the maximal channel constriction $h_{m}^{(s)}=\min _{x}(h^{(s)})$ with corresponding $x$ position $x_{m}^{(s)}$ .

2.2 Linear stability global eigensolver

We consider the stability of a static solution of (2.9) by considering the growth of a small perturbation to this state, writing

where $\unicode[STIX]{x1D703}\ll 1$ and $\tilde{h}(x;\unicode[STIX]{x1D70E})$ and $\tilde{q}(x;\unicode[STIX]{x1D70E})$ are functions of $x$ and $^{\ast }$ denotes a complex conjugate. Substituting into the full governing equations (2.8) and neglecting terms of $O(\unicode[STIX]{x1D703}^{2})$ we obtain the linearised equations

subject to linearised boundary conditions

Equations (2.11) constitute a global eigenvalue problem for $\unicode[STIX]{x1D70E}$ and the corresponding eigenfunctions $\tilde{q}(x;\unicode[STIX]{x1D70E})$ and $\tilde{h}(x;\unicode[STIX]{x1D70E})$ , which must be solved numerically. To form a solution to the eigenvalue problem (2.11) we construct a global eigensolver by discretising the collapsible segment of the channel $x$ onto a uniform mesh $x_{j}=j\unicode[STIX]{x0394}x$ ( $j=0,\ldots ,N-1$ ), where $\unicode[STIX]{x0394}x=1/(N-1)$ . We then discretise the perturbation eigenfunctions onto this mesh by writing $\tilde{h}_{j}=\tilde{h}(x_{j};\unicode[STIX]{x1D70E})$ and $\tilde{q}_{j}=\tilde{q}(x_{j};\unicode[STIX]{x1D70E})$ ( $j=0,\ldots ,N-1$ ) and form the vector of unknown coefficients $\boldsymbol{a}=[\tilde{h}_{0},\ldots ,\tilde{h}_{N-1},\tilde{q}_{0},\ldots ,\tilde{q}_{N-1}]^{\text{T}}$ . By expressing spatial derivatives using second-order centred finite difference approximations, we hence express (2.11) as an eigenvalue problem in the form

where $\unicode[STIX]{x1D648}$ is a matrix of size $2N\times 2N$ . The eigenvalues ( $\unicode[STIX]{x1D70E}$ ) and corresponding eigenvectors ( $\boldsymbol{a}$ ) can be calculated numerically using Matlab in a straightforward manner. For computation of the energy budget below we normalise the eigenfunction so that $\max (\text{Re}(\tilde{q}(x;\unicode[STIX]{x1D70E})))=1$ unless specified otherwise. Results of the linear stability analysis are described in §§ 4.1–4.2. In simulations we use up to a maximum of $N=1600$ grid points to resolve the highly constricted static shape and associated eigenfunctions. Simulations for larger $N$ became impractical and were not deemed necessary. This global eigensolver is similar in spirit to that constructed by Luo et al. (Reference Luo, Cai, Li and Pedley2008) (see also Liu et al. (Reference Liu, Luo, Cai and Pedley2009, Reference Liu, Luo and Cai2012)) using a finite element method and an Arnoldi algorithm for computing eigenvalues, but their two-dimensional eigensolver requires extremely large matrices and is impractical for large-scale parameter sweeps. Since the present approach is one-dimensional it has the advantage of much reduced computational cost permitting a unified overview of the parameter space.

2.3 Fully nonlinear simulations

Nonlinear, time-dependent simulations of (2.8) can be constructed using the semi-implicit finite difference method described and used previously by Stewart et al. (Reference Stewart, Waters and Jensen2009). In particular, we discretise spatial derivatives using fourth-order finite differences and time derivatives using first-order finite differences, representing the highest spatial derivative in each equation implicitly with all other terms treated explicitly. To capture the highly complicated ‘slamming’ oscillations described below we use an increased temporal and spatial grid in regions where the membrane gets drawn close to the rigid wall. We have also thoroughly checked for convergence of the scheme using a variety of spatial and temporal grids up to a maximum of 3000 equally spaced spatial grid points. Results of the fully nonlinear simulations are described in § 7. We typically characterise the numerical simulations by considering the temporal evolution of the maximal channel constriction $h_{m}(t)=\min _{x}(h(x,t))$ .

2.4 Energy budget

Following Stewart et al. (Reference Stewart, Waters and Jensen2009), we compute the corresponding energy budget for the flow by multiplying (2.8) by $6q/(5h)$ and integrating over the channel length and width. This energy budget can be written in the form

where

where ${\mathcal{K}}$ measures the rate of change of kinetic energy, ${\mathcal{F}}$ is the net kinetic energy flux extracted from the mean flow, ${\mathcal{P}}$ is the rate of working of upstream pressure forces, ${\mathcal{E}}$ is the rate of working of pressure forces across the membrane, ${\mathcal{D}}_{S}$ is the rate of working of viscous forces and ${\mathcal{D}}_{A}$ is the additional kinetic energy flux induced by the flow profile assumption. The terms in the fully nonlinear energy budget for periodic oscillations are computed in § 7.3. It should be noted that ${\mathcal{P}}$ is proportional to the upstream pressure $p_{u}$ , which can in turn be expressed in terms of the external pressure on the membrane and the membrane curvature at the upstream end of the flexible segment.

2.5 Energy budget for static solutions

For a static membrane (described in § 2.1) the only non-trivial contributions to the energy budget can be expressed as (steady components denoted with a superscript $(s)$ ),

where

Using the downstream boundary condition (2.9c ) this static energy balance can be rewritten as

this relation is employed to close the reduced system of equations constructed in the limit of large $p_{e}$ in § 6.

2.6 Rate of working of viscous dissipation

To elucidate the energy budget of neutrally stable oscillations, we consider the rate of working of viscous forces (2.13f ) in more detail. In particular, we consider the time average over a period of oscillation, denoted by an overbar. We consider the perturbation expansion posed in (2.10) continued to the next order in $\unicode[STIX]{x1D703}$ in the form

Using this expansion we approximate ${\mathcal{D}}_{S}$ (2.13f ) by expanding in powers of $\unicode[STIX]{x1D703}$ . At leading order $\overline{{\mathcal{D}}_{S}}={\mathcal{D}}_{S}^{(s)}$ . At the following order, the variables have components which are proportional to either $\exp (\unicode[STIX]{x1D70E}t)$ or $\exp (\unicode[STIX]{x1D70E}^{\ast }t)$ , so taking the time average over a period of neutrally stable oscillation ( $\text{Re}(\unicode[STIX]{x1D70E})=0$ ) implies that the $O(\unicode[STIX]{x1D703})$ contribution to $\overline{{\mathcal{D}}_{S}}$ must be zero. Conversely, at $O(\unicode[STIX]{x1D703}^{2})$ , the checked variables consist of time-independent terms (denoted as $\overline{\check{h}}$ and $\overline{\check{q}}$ for $h$ and $q$ , respectively) as well as time-dependent components which are proportional to either $\exp (2\unicode[STIX]{x1D70E}t)$ or $\exp (2\unicode[STIX]{x1D70E}^{\ast }t)$ . Hence, when averaging over a period of neutrally stable oscillation, we obtain a non-trivial contribution to the energy budget at $O(\unicode[STIX]{x1D703}^{2})$ . A similar approach was previously considered by Stewart et al. (Reference Stewart, Heil, Waters and Jensen2010) to compute the perturbation energy budget of ‘sloshing’ oscillations driven by fixed upstream pressure.

Hence, we expand the time-averaged rate of working of viscous dissipation in the form

where

consists of terms which involve the eigenfunctions $\tilde{h}$ and $\tilde{q}$ (and their complex conjugates) determined above, while

consists only of time-independent terms arising at second order in $\unicode[STIX]{x1D703}$ . We note that $\overline{\check{q}}$ must be zero to satisfy conservation of mass over a period of neutrally stable oscillation, so $\overline{\check{{\mathcal{D}}}_{S}}$ depends on $\overline{\check{h}}$ alone. The quantity $\overline{\tilde{{\mathcal{D}}}}_{S}$ can be easily computed from eigenfunctions obtained from the global stability eigensolver discussed in § 2.2. However, solution of the $O(\unicode[STIX]{x1D703}^{2})$ perturbation problem would be required to compute $\overline{\check{{\mathcal{D}}}_{S}}$ , which is not considered here. We use this decomposition of $\overline{{\mathcal{D}}_{S}}$ in § 7.3 to elucidate the mechanism of energy transfer sustaining the self-excited oscillations.

4.1 Typical eigenvalue spectra

Having characterised the static membrane profiles, we now use the global eigensolver described in § 2.2 to consider their stability to time-dependent perturbations. Two typical eigenvalue spectra are plotted in figure 3(a), plotting $\text{Re}(\unicode[STIX]{x1D70E})$ against $\text{Im}(\unicode[STIX]{x1D70E})$ for the static state plotted in the inset. Eigenvalues appear as complex conjugate pairs, which are all stable ( $\text{Re}(\unicode[STIX]{x1D70E})<0$ ) for sufficiently small Reynolds numbers. As the Reynolds number increases a pair of eigenvalues will cross the $\text{Re}(\unicode[STIX]{x1D70E})$ axis with non-zero imaginary part, indicative of a Hopf bifurcation. Two cases are illustrated for $T=2$ , $p_{e}=10$ and $L_{2}=3$ in figure 3(a): $R=61$ , where all eigenvalues are stable (filled black circles), and $R=65$ , where two conjugate pairs of eigenvalues are unstable (crosses). In figure 3(b–d) we plot the corresponding eigenfunctions $\tilde{h}(x;\unicode[STIX]{x1D70E})$ of the three most unstable modes for $R=65$ ; the structure of these eigenfunctions is explained in more detail below.

Figure 3.

(a) Eigenvalue spectra for $T=2$ , $p_{e}=10$ and $L_{2}=3$ , with $R=61$ (filled circles) and $R=65$ (crosses). The pairs of modes are labelled (i)–(vi) in order of increasing stability. Inset shows the corresponding static shape of the wall for $R=61$ (dashed) and $R=65$ (solid). Corresponding eigenfunctions $\tilde{h}(x;\unicode[STIX]{x1D70E})$ of the three most unstable modes for $R=65$ : (b) mode (i); (c) mode (ii); (d) mode (iii). Eigenfunctions have been normalised so that $\max (\text{Re}(\tilde{q}))=1$ . The dotted line in (b–d) denotes the position $x_{m}^{(s)}$ .

4.2 Neutrally stable oscillations

Neutral stability curves (along which $\text{Re}(\unicode[STIX]{x1D70E})=0$ ) in the parameter space spanned by the membrane tension and critical Reynolds number are shown in figure 4(a) for $p_{e}=10$ and $L_{2}=3$ . In particular, we illustrate the neutral stability curves for the five most unstable modes identified in the spectrum in figure 3. These modes are labelled (i) to (v). We have verified that for these parameter values the primary global instability of the system is always either to the mode labelled (i) or (iia). It should be noted that the mode labelled (ii) becomes unstable in two separate regions across the parameter space, labelled (iia) and (iib), respectively. Interestingly, the neutral curve for mode (iib) exhibits a low-frequency oscillation for low tensions (consistent with the framework described in § 5) but a high-frequency oscillation for large tensions (consistent with the framework described in § 6). Although these neutral curves intersect one another, they do so with distinct frequencies, so there is no mode interaction. In addition, in figure 4(b) we illustrate regions of parameter space where the system is unstable to oscillations using hatching, as well as highlighting the corresponding regions of parameter space where multiple steady states exist (dark shaded region). The corresponding oscillation frequency of each mode is shown in figure 4(c); for each mode the oscillation frequency increases with membrane tension, which is exploited in § 6.

Previous studies of collapsible channel flow (e.g. Luo et al. Reference Luo, Cai, Li and Pedley2008; Stewart et al. Reference Stewart, Waters and Jensen2009) characterise these modes by the number of turning points in the eigenfunction of the membrane shape, $\tilde{h}(x;\unicode[STIX]{x1D70E})$ . Hence, according to this convention for $T=2$ the mode labelled (i) would be mode-4, the mode labelled (iia) would be mode-5 and the mode labelled (iii) would be mode-6. However, for lower tensions the mode labelled (i) has two extrema (see figure 7 a) and is thus mode-2, so the shape of the eigenfunction changes along the neutral curve. In this study it emerges that for large tensions a more convenient labelling convention is to characterise the eigenfunction relative to the point of greatest channel constriction, $x_{m}^{(s)}$ (position denoted with a dotted line on the panels in figure 3 b–d). We characterise the eigenfunctions of neutrally stable modes in more detail using the asymptotic model of § 6.

Figure 4.

Neutrally stable modes for $p_{e}=10$ and $L_{2}=3$ as a function of the membrane tension $T$ : (a) the critical value of the Reynolds number for instability of the six most unstable modes labelled (i) to (v); (b) regions of the $(T,R)$ parameter space where the system has multiple coexisting static states (dark shaded region) and regions where the highly collapsed basic state is unstable to oscillations (hatching); (c) corresponding frequency of the neutrally stable modes. The open circles in (a) correspond to limit (saddle) points associated with a change in the number of static configurations, while the crosses represent the leading-order prediction for the critical Reynolds number from the low-frequency theory in § 5. The labels $(m,n)$ follow from the asymptotic description of the eigenfunction of the wall shape in the limit of large $T$ described in § 6, where $m$ is the number of turning points upstream of the constriction and $n$ is the number of turning points downstream.

It is also instructive to consider the onset of instability in terms of the corresponding static configuration of the membrane. In figure 2(a) we indicate with a filled circle the critical value of the membrane tension below which each branch of static solutions becomes unstable to self-excited oscillations. For this value of the external pressure, configurations approaching the uniform state (along the upper branch as $T\rightarrow \infty$ ) are always stable for large tensions. This is consistent with the study of Xu et al. (Reference Xu, Billingham and Jensen2013), who found that the uniform static configuration (maintained with an external pressure gradient) is always stable beyond a critical value of the membrane tension $T\approx 0.121$ .

Of most pertinent interest in the system is the first onset of self-excited oscillations, termed the primary oscillatory instability. For increasing $T$ , the system exhibits three distinct behaviours. Each of these is now discussed in turn.

6.1 Region 1

To compute an approximation to the static membrane shape in region 1 ( $0\leqslant x\leqslant x_{0}$ , denoted with the subscript 1), we expand the static membrane profile according to

In region 1 the leading-order static governing equations (2.9) become

subject to boundary conditions,

This system can be solved to obtain

where $\overline{A}_{1}$ is a constant which is determined by matching into region 3.

Similarly, to construct an approximation to the perturbation eigenfunction across region 1 we expand the dependent variables according to

To leading order across region 1 the linearised governing equations (2.11) can be written entirely in terms of $\hat{q}_{1}$ in the form

subject to boundary conditions

This represents a leading-order balance between the unsteady inertial forces and the restoring force of membrane tension. It should be noted that only $\hat{\unicode[STIX]{x1D70E}}^{2}$ appears in the governing equations. Hence, at neutral stability ( $\hat{\unicode[STIX]{x1D70E}}=\text{i}\hat{\unicode[STIX]{x1D6FD}}$ ) the equations for real and imaginary parts of $\hat{q}$ decouple entirely across this region.

We solve the boundary value problem (6.5a ) using a shooting method based on the ode15s solver implemented in Matlab; simulations are initiated from $x=x_{0}$ , using the two specified boundary conditions (6.5c ). In addition, we adjust the values of the boundary conditions $\hat{q}_{1,xx}(x_{0}^{-})$ and $\hat{q}_{1,xxx}(x_{0}^{-})$ to isolate solutions which obey the appropriate boundary conditions at $x=0$ (6.5b ). However, the boundary value problem (6.5a ) is locally singular at $x=x_{0}$ since ${\hat{h}}_{1}^{(s)}(x_{0})=0$ . To initiate a numerical simulation we shoot from a small distance $\unicode[STIX]{x1D6FF}$ upstream of $x_{0}$ ( $\unicode[STIX]{x1D6FF}\ll 1$ ) towards the upstream boundary $(x=0)$ . We approximate the value of $\hat{q}_{1}$ and its derivatives at this point using a Taylor series expansion about $x_{0}$ . In simulations we found that a six-term Taylor series with $\unicode[STIX]{x1D6FF}=10^{-3}$ was sufficient to ensure convergence.

The corresponding perturbation membrane profile across region 1 is given by ${\hat{h}}_{1}=-\hat{q}_{1,x}/\hat{\unicode[STIX]{x1D70E}}$ . At neutral stability, assuming $\hat{q}_{1}$ is purely real implies that ${\hat{h}}_{1}$ will be purely imaginary.

6.2 Region 2

In a similar manner to region 1, across region 2 ( $x_{0}\leqslant x\leqslant 1$ , denoted with the subscript 2), we expand the static membrane profile according to

In region 2 the leading-order static governing equations (2.9) become

subject to boundary conditions

This system can be easily solved to obtain

To construct an approximation to the perturbation eigenfunction across region 2 we expand the dependent variables according to

The perturbation governing equations (2.11) in region 2 take the leading-order form

subject to

This represents a similar balance of terms to region 1.

Again this problem is locally singular at $x=x_{0}$ and we employ an analogous numerical method to that described in § 6.1 to construct a solution; we shoot from a small distance $\unicode[STIX]{x1D6FF}$ downstream of the constriction $x_{0}$ using the boundary conditions (6.9b ) towards the downstream boundary at $x=1$ , adjusting the value of $\hat{q}_{2,xx}(x_{0}^{+})$ and $\hat{q}_{2,xxx}(x_{0}^{+})$ to identify solutions which satisfy the boundary conditions (6.9c ).

The corresponding perturbation membrane profile across region 2 is then given by ${\hat{h}}_{2}=-\hat{q}_{2,x}/\hat{\unicode[STIX]{x1D70E}}$ . As before, at neutral stability we assume $\hat{q}_{2}$ is purely real, so ${\hat{h}}_{2}$ must be purely imaginary.

6.3 Region 3

In region 3, close to the point of maximal constriction at $x=x_{0}$ , we rescale the variables according to

Hence, across region 3, the leading-order governing equation for the static membrane shape (2.9) becomes

representing an inviscid balance between convective inertia and the restoring force of membrane tension. This equation is subject to continuity of slope in the far field with regions 1 and 2, in the form

Following Xu & Jensen (Reference Xu and Jensen2015), we integrate (6.11a ) twice to form an explicit expression for the membrane slope in the form

where ${\hat{H}}_{3}^{(s)}(0)=H_{0}$ is an unknown parameter. Following Xu & Jensen (Reference Xu and Jensen2015), we can integrate once more to obtain an implicit expression for the membrane shape across region 3,

For use below, following Xu & Jensen (Reference Xu and Jensen2015) it can also be shown that

Furthermore, to construct an approximation to the perturbation eigenfunction across region 3 we expand the dependent variables according to

The leading-order perturbation equations across region 3 become

subject to appropriate matching conditions into regions 1 and 2. However, noting that these equations are symmetric in $X$ we find below that we do not need to solve the perturbation equations in this region explicitly to construct the outer approximation to the eigenfunctions.

6.4 The composite static solution

Matching the expression for ${\hat{H}}_{3,X}^{(s)}$ as $X\rightarrow -\infty$ into region 1 we obtain the constant $\overline{A}_{1}$ in the form

while matching the expression for ${\hat{H}}_{3,X}^{(s)}$ as $X\rightarrow \infty$ into region 2 we obtain a constraint for $H_{0}$ as a function of $x_{0}$ ,

This fully specifies the leading-order static solution across all three regions for a given $x_{0}$ . Using the value of $x_{0}$ calculated in § 6.5 and the critical value of Reynolds number $\hat{R}$ calculated in § 6.6, we thus plot the approximation to the membrane shape across the three regions in figure 7 for three examples (figure 7 d,g,j), showing excellent qualitative agreement with full numerical simulations in the first two cases; a discussion of the discrepancy for mode (2,1) is given in § 6.7.

6.5 The oscillation frequency and point of maximal constriction

Since we are only concerned with neutrally stable oscillations, we express $\hat{\unicode[STIX]{x1D70E}}^{2}=-\hat{\unicode[STIX]{x1D6FD}}^{2}$ and conduct all calculations in terms of $\hat{\unicode[STIX]{x1D6FD}}$ . Given that $\hat{\unicode[STIX]{x1D6FD}}^{2}$ is purely real at neutral stability, the form of the governing equations (6.5a ), (6.9a ) mean that the real and imaginary parts of $\hat{q}$ can be considered independently at this order.

Across region 3 we observe that the membrane profile is symmetric to leading order, so that the limiting slopes from regions 1 and 2 are equal and opposite,

Since the solution to the linearised system (2.11) can be normalised arbitrarily, we set $\hat{q}_{xx}(x_{0}^{+})=1=-\hat{q}_{xx}(x_{0}^{-})$ without loss of generality (since $\hat{q}_{xx}$ is proportional to ${\hat{h}}_{x}$ through (2.11a )). For a given value of $x_{0}$ we compute the corresponding values of $\hat{\unicode[STIX]{x1D6FD}}$ and $\hat{q}_{xxx}(x_{0}^{-})$ ( $\hat{q}_{xxx}(x_{0}^{+})$ ) such that the solution satisfies the upstream (downstream) boundary conditions, from which we trace out curves in the parameter space spanned by $x_{0}$ and $\hat{\unicode[STIX]{x1D6FD}}$ ; each curve represents a different number of extrema in the profile of $\tilde{h}$ denoting $m$ ( $n$ ) turning points upstream (downstream) of the point of maximal constriction. In figure 8, typical curves for $m=2$ , $m=3$ and $m=4$ are shown as solid lines while curves for $n=1$ , $n=2$ and $n=3$ are shown as dashed lines, computed for $\hat{T}=0.2$ and $L_{2}=3$ . Numerical tests indicate there is no curve for $m=1$ . Points of intersection between upstream (solid) and downstream (dashed) curves in figure 8 represent global roots of the full system. These global crossing points correspond to distinct modes with a particular $(m,n)$ ; the critical values of $x_{0}$ and $\hat{\unicode[STIX]{x1D6FD}}$ for the first nine modes are listed in table 1, along with the corresponding shooting parameters.

Figure 8.

Parameter space spanned by the position of the membrane minimum $x_{0}$ against the corresponding oscillation frequency $\hat{\unicode[STIX]{x1D6FD}}$ for $\hat{T}=0.2$ and $L_{2}=3$ . Solid (dashed) lines represent solutions which satisfy the upstream (downstream) boundary conditions. Crossing points represent global roots of the asymptotic system.

Table 1.

Shooting parameters for the nine eigenmodes shown on figure 8 computed for $\hat{T}=0.2$ and $L_{2}=3$ .

The asymptotic prediction of the critical frequencies for modes (3,2), (4,2) and (2,1) are plotted against the full numerical solutions as filled circles in figures 4(c) and 10(b). We find excellent agreement with the computed frequencies, with the mode (2,1) corresponding to neutral curve (iib), the mode (3,2) corresponding to neutral curve (iv) and the mode (4,2) corresponding to the neutral curve (v). In addition, the corresponding asymptotic prediction of the maximal constriction ( $p_{e}^{-1}H_{0}$ ) and the point of greatest constriction ( $x_{0}$ ) are shown in figures 10(c) and 10(d), respectively, against the full numerical computations. In both cases the curves (3,2) and (4,2) show excellent agreement with the numerical computations for neutral curves (iv) and (v). The curve (2,1) shows the correct trend as neutral curve (iib), but is slightly separated. The origin of this discrepancy is explored in § 6.7.

6.6 The critical Reynolds number

Applying the scalings for regions 1 to 3 into the energy budget for static solutions (2.17), we obtain the leading-order relationship

so the rate of working of pressure forces along the compliant segment of the channel is balanced by the viscous losses close to the point of greatest channel constriction. Substituting the analytical expressions for $h_{1}^{(s)}$ , $h_{2}^{(s)}$ constructed above (6.3d ) and (6.7e ) and the integral (6.14), we obtain the constraint

For a given $x_{0}$ and $\hat{\unicode[STIX]{x1D6FD}}$ computed using the method outlined above, the corresponding critical Reynolds number can be computed using (6.21). These predictions are listed in table 1 and the corresponding predictions of the neutral curves are shown in figures 4(a) and 10(a) for modes (2,1), (3,2) and (4,2). As $p_{e}\rightarrow \infty$ the neutral curves each show the correct trend and approach the numerically computed curves, although the values of $p_{e}$ considered here are not sufficiently large to demonstrate a strong overlap (as discussed above, we found that an unfeasible number of grid points were needed to resolve the static profile for larger values of $p_{e}$ ). There is a more marked discrepancy for the mode labelled (2,1), which is discussed below. For some of the roots identified in table 1 the corresponding critical Reynolds number is negative, indicating an unphysical solution for those $(m,n)$ and that choice of parameters $\hat{T}$ and $L_{2}$ .

6.7 Outer approximation to the eigenfunctions

Given the approximate eigenfunction profiles constructed across regions 1 and 2 in § 6.5, the outer asymptotic eigenfunctions for three particular examples are shown in figure 7(d–l) for $\hat{T}=0.2$ and $L_{2}=3$ . Note that the outer asymptotic eigenfunction has been normalised in the same way as the numerically computed eigenfunctions for ease of comparison. In particular, we consider mode (3,2) in figure 7(e,f), where $R_{c}\approx 338.44$ ( $p_{e}\hat{R}_{c}\approx 216.30$ ), mode (4,2) in figure 7(h,i), where $R_{c}\approx 677.65$ ( $p_{e}\hat{R}_{c}\approx 532.15$ ) and mode (2,1) in figure 7(k,l), where $R_{c}\approx 320.09$ ( $p_{e}\hat{R}_{c}\approx 605.69$ ). For modes (3,2) and (4,2) we observe excellent qualitative agreement between the numerically computed eigenfunctions and the outer asymptotic approximation. However, for mode (2,1) the agreement is weaker, exacerbated by the large discrepancy in the estimate of the point of maximal constriction (figure 10 d). A possible explanation for this discrepancy lies in the shape of the eigenfunction downstream of the constriction. In this region, unlike for the other examples shown, the real and imaginary parts of $\tilde{q}$ are comparable (figure 10 d,ii), so our assumption that $\tilde{q}$ is purely real is not appropriate. Furthermore, the slopes of the perturbation membrane shape upstream and downstream of the constriction are clearly not equal and opposite (figure 10 d,iii), in contradiction of our leading-order asymptotic theory. To improve the prediction of this mode would require coupling of the real and imaginary parts of the eigenfunction at leading order, which is not considered here.

7.1 Primary global instability

We use nonlinear simulations of (2.8) to trace the behaviour along two lines in parameter space with $p_{e}=10$ and $L_{2}=3$ , the first for $T=2$ , where the collapsed basic state exists for all $R$ (§ 7.1), and the second for $T=2.7$ , where the collapsed basic state (§ 7.2) does not exist for sufficiently small $R$ .

Figure 9.

Fully nonlinear simulations for $p_{e}=10$ and $L_{2}=3$ : (a) bifurcation diagram of $h_{m}$ against $R$ for $T=2$ , tracing the maximum (filled circles) and minimum (crosses) of $h_{m}(t)$ over a period of oscillation; time traces of the membrane minimum $h_{m}$ for (b) $T=2$ , $R=61$ , where the response is asymptotically stable; (c) $T=2$ , $R=65$ , where the oscillatory response has undergone a period doubling bifurcation; (d) $T=2.7$ , $R=120$ where the oscillation again saturates the almost uniform static state; (e) $T=2.7$ , $R=130$ where the oscillation saturates the almost uniform static state; (f) $T=2.7$ , $R=170$ where the oscillation saturates to a large-amplitude limit cycle with ‘slamming’. The triangles in (a) illustrate the mean position of the membrane minimum over a period, while the thick black lines indicate the range of values for the maximum and minimum of $h_{m}$ for the chaotic oscillations.

We examine the primary global instability close to the neutral point for $T=2$ , $p_{e}=10$ for $L_{2}=3$ , where the collapsed basic state becomes unstable to an oscillatory instability (via a Hopf bifurcation) as the Reynolds number increases through the threshold $R=R_{c}\approx 62.24$ . In simulations we employ an initial condition

where $A=0.01$ in simulations. In figure 9(a) we numerically trace a bifurcation diagram by plotting the minimal and maximal values of $h_{m}$ over a period of fully developed oscillation as a function of $R$ . Figure 9(a) indicates that the Hopf bifurcation is supercritical. For $R<R_{c}$ the static membrane is asymptotically stable (figure 9 b). However, as the Reynolds number increases the oscillation undergoes a period doubling bifurcation (evident for $R=65$ , see figure 9 c) and eventually becomes chaotic (evident for $R=66$ , $R=67$ , where the thick black lines indicate the range of values for the maximum and minimum of $h_{m}$ ). Also shown on figure 9(a) is the mean position of the membrane minimum over a period of oscillation, denoted $\overline{h_{m}}$ (illustrated as open triangles connected by a dotted line) compared to the minimal value of the corresponding static membrane (solid line). We observe that the mean minimal channel width is increased slightly compared to the static solution, which has consequences for the energy transfer mechanism sustaining the oscillation (see § 7.3).

7.2 Fully nonlinear ‘slamming’ oscillations

We examine the secondary global instabilities by considering the fully nonlinear behaviour close to a region where the collapsed basic state does not exist, in particular for $T=2.7$ , $p_{e}=10$ for $L_{2}=3$ . For $R>R_{b}\approx 127.8$ we use the initial condition (7.1a,b ). For lower Reynolds numbers ( $R<R_{b}$ ) the highly collapsed steady state does not exist. Starting close to the almost uniform static state the system is asymptotically stable (not shown). Hence, we choose to start the simulations close to the highly collapsed static state computed for some $R>R_{b}$ (here we choose an initial state generated for $R=130$ ). For $R=120$ the system oscillates initially but eventually saturates towards the almost uniform static state, as in figure 9(d). For $R>R_{b}$ the perturbation growth becomes increasingly oscillatory and irregular as $R$ increases, where there are two possible attractors for long times. For some parameter choices inside the dark shaded region in figure 4(b) (where there are two concurrent static states) the system oscillates initially but eventually saturates to the almost uniform static profile, as shown in figure 9(e) for $R=130$ ; for other parameter choices the system exhibits a fully developed oscillatory (typically chaotic) limit cycle, where periodically the membrane is drawn very close to the rigid wall before recovering (see example for $R=170$ in figure 9 f). These so-called ‘slamming’ oscillations have been observed previously (Stewart Reference Stewart2010; Stewart et al. Reference Stewart, Heil, Waters and Jensen2010; Xu et al. Reference Xu, Billingham and Jensen2013) and can be deconstructed using matched asymptotic expansions (Xu & Jensen Reference Xu and Jensen2015). As the Reynolds number increases the ‘slamming’ oscillation eventually becomes dominant. It is interesting to note that in this case linear instability to oscillations does not necessarily result in an unsteady outcome for long times, similar to predictions of three-dimensional flow in a collapsible tube by Heil & Boyle (Reference Heil and Boyle2010).

7.3 Energy budget of the primary global instability

Finally, we examine the energy budget of the self-excited oscillations to elucidate the driving mechanism. We define the time average of a quantity $f(x,t)$ over a period of oscillation $T_{f}$ according to

In particular we consider a parameter sweep along the line $T=2$ for $p_{e}=10$ , $L_{2}=3$ where the critical Reynolds number $R_{c}\approx 62.24$ corresponds to the mode labelled (i) in figure 4. The terms in the full energy budget (2.13a ) are computed for various $R>R_{c}$ from the nonlinear simulations and their time average calculated over a period of fully saturated oscillation and listed in table 2. Also listed in the table are the terms in the energy budget for the corresponding static state, denoted with the superscript $(s)$ .

We compute the energy budget of the oscillation by subtracting the energy corresponding to the static state, computing the excess work done by upstream pressure forces over a period ( $\overline{{\mathcal{P}}^{(e)}}=\overline{{\mathcal{P}}}-{\mathcal{P}}^{(s)}$ ) and the excess rate of working of viscous dissipation ( $\overline{{\mathcal{D}}^{(e)}}=\overline{{\mathcal{D}}_{S}}-{\mathcal{D}}^{(s)}$ ). It emerges that over a period the dominant energy balance is between these excess energies, $\overline{{\mathcal{P}}^{(e)}}$ and $\overline{{\mathcal{D}}^{(e)}}$ , while the net kinetic energy flux extracted from the mean flow ( $\overline{{\mathcal{F}}}$ ) and the kinetic energy flux induced by the flow profile assumption ( $\overline{{\mathcal{D}}_{A}}$ ) are both very small in comparison (see table 2). This is in direct agreement with the observations of Liu et al. (Reference Liu, Luo and Cai2012) in a flow-driven collapsible channel system.

As discussed in § 2.4, we further decompose $\overline{{\mathcal{D}}_{S}}$ using the global linear stability eigensolver, where (2.19) implies

for $\unicode[STIX]{x1D703}\ll 1$ . In the fully nonlinear computations we find $\overline{{\mathcal{D}}^{(e)}}<0$ for all cases tested (see table 2). However, computing $\overline{\hat{{\mathcal{D}}}_{S}}$ for these parameter values at criticality ( $T=2$ , $p_{e}=10$ , $L_{2}=3$ , $R_{c}\approx 62.24$ ) using the perturbation eigenfunctions $\tilde{h}$ and $\tilde{q}$ we find that $\overline{\hat{{\mathcal{D}}}_{S}}>0$ . Hence, dissipative forces in the perturbation are consuming energy, so the perturbation dissipation in the mean flow, $\overline{\breve{{\mathcal{D}}}_{S}}$ , must be negative. Hence, the presence of the oscillation reduces the overall dissipation of the mean flow (through the work of nonlinear Reynolds stresses) and this excess energy drives the oscillation. However, this reduction of the dissipation is not achieved by modifying the mean flow along the channel as we need $\overline{q}(1,t)=1$ (i.e.  $\overline{\check{q}}=0$ ) to satisfy conservation of mass over a period of neutrally stable oscillation. Instead this reduction in dissipation is achieved by decreasing the mean maximal constriction of the channel over a period of neutrally stable oscillation, so $\overline{h_{m}}>h_{m}^{(s)}$ . This is illustrated in the fully nonlinear bifurcation diagram in figure 9(a) and the computed values of $\overline{h_{m}}$ and $h_{m}^{(s)}$ are listed in table 2.

A similar mechanism of energy transfer is operating for the mode labelled (iia), corresponding to the most unstable member of the family of high-frequency oscillations explored in § 6.

Table 2.

Numerical values of time-averaged quantities in the fully nonlinear energy equation (2.13) and corresponding static energy equation (2.14) for simulations with fixed $T=2$ , $p_{e}=10$ and $L_{2}=3$ .