## Abstract

The effects of wall inertia on instabilities in a collapsible channel with a long finite-length flexible wall containing a high Reynolds number flow of incompressible fluid are studied. Using the ideas of interactive boundary layer theory, the system is described by a one-dimensional model that couples inviscid flow outside the boundary layers formed on the channel walls with the deformation of the flexible wall. The observed instability is a form of flutter, which is superposed on the behaviour of the system when the wall mass is neglected. We show that the flutter has a positive growth rate because the fluid loading acts as a negative damping in the system. We discuss these findings in relation to other work on self-excited oscillations in collapsible channels.

## 1. Introduction

Self-excited oscillations in a collapsible channel have been studied widely in the fluid-mechanical literature, because of their relevance to various physiological phenomena, such as wheezing in lung airways and Korotkov sounds in blood-pressure measurement. Similar instabilities are readily observed in experimental collapsible-tube systems, such as the Starling resistor [1], but are yet to be fully understood in all flow regimes [2]. Moreover, although collapsible channels are inherently global systems with self-excited oscillations often being driven or considerably modified by the presence of finite boundaries [3], they are undoubtedly linked to the local instabilities in homogeneous domains with flexible walls [4–7].

In the simplest collapsible channel model, a finite-Reynolds-number fluid flow is driven through a two-dimensional channel that has a finite segment of one of the walls replaced with a prestretched elastic membrane under external pressure. Depending on the system parameters and the controlled boundary conditions, different types of oscillations that correspond to different modes, or number of local extrema in the wall shape, are observed. When the pressure drop is controlled in the system, it typically loses stability through mode-1 ‘sloshing’ oscillations, described by Jensen & Heil [8] and Stewart *et al.* [9]. The mechanism for sustaining these oscillations relies on the kinetic energy flux into the system being greater than the energy loss due to viscous dissipation. Hence ‘sloshing’ is observed only in channels whose downstream rigid segments are longer than the upstream rigid segments. At lower values of the wall tension, these oscillations are replaced by large-amplitude ‘slamming’ oscillations, during which the membrane briefly almost touches the opposite wall [9]. On the other hand, when the flow rate is the controlled boundary condition, the self-excited oscillations appear as a consequence of instability of the steady solution due to small-amplitude mode-2 perturbations [10,11]. If the flexible wall of the collapsible channel is long, but similar in length to the downstream rigid section, then an oscillatory mode arises through an interaction between two static eigenmodes Xu *et al.* [12]. Crucially however neither of the classical scenarios includes significant wall inertia, which is believed to complicate the dynamical picture even further. For example, for a non-zero wall mass Luo & Pedley [13] observed a higher frequency flutter instability superimposed on the original large-amplitude, low-frequency oscillations. Unlike the oscillations of the massless wall, the flutter grew with time from a small amplitude until it completely dominated the original slower mode. Therefore, one anticipates a whole new range of mechanisms operating when the structural inertia is accounted for as well.

One such mechanism was explored by Huang [14], who studied a finite-length channel conveying a Poiseuille flow with a segment of one of the walls replaced by a heavy tensioned membrane. By focusing on the linear stability of the uniform Poiseuille flow, the author identified eigenmodes related to flutter, which were stabilized by the structural damping in the system, but destabilized by the fluid viscosity. The asymmetric loading created by the fluid on the membrane led to oscillations with different preferential modes in the membrane shape depending on the controlled boundary conditions.

Another flutter mechanism in collapsible channels was studied by Mandre & Mahadevan [15], who explained it as a 1:1 resonance involving a coincidence of two frequencies in a dynamical system. However, this type of flutter was limited to systems dominated by the normal mode oscillations of elastic components, but also susceptible to weak fluid forcing with time scale much smaller than the natural time scale of the elastic system oscillations.

Here, we present a different form of flutter in a collapsible channel, which relies on coupling between the inviscid flow outside the boundary layers formed on the channel walls and the deformation of the flexible heavy membrane under longitudinal tension. We study this model in the context of the previous works by Luo & Pedley [10,13]. We start by presenting our model based on interactive boundary layer theory in §2. In §3 we present our results, and discuss them in §4.

## 2. The model

We focus on a high Reynolds number flow (*Re*≫1) of incompressible viscous fluid, with mean velocity *U*_{0}, density *ρ* and kinematic viscosity *ν*, in a two-dimensional channel of width *a* with a long slender indentation in the lower wall of a finite length λ*a*, λ≫1, and amplitude *aεF*, *ε*≪1 (figure 1). The flow rate at the entrance is taken to be fixed. Under the following non-dimensionalization:
*Re*^{1/7}≤λ≪*Re* and *Re*^{−2/7}≪*ε*≪1 [16]. Here, *u*_{0}(*y*)=6*y*(1−*y*) is the Poiseuille velocity profile, −12λ*x*/*Re* is a small term required to drive the Poiseuille flow, −*εA*(*x*,*t*) is the lateral displacement of streamlines in the core flow, which would in the absence of the indentation stay straight, and *P* is the unknown contribution to the pressure *p*. From (2.3), it can also be seen that
*σ*=6/5λ^{2}*ε* represents the importance or otherwise of the cross-stream pressure gradient. Similar systems have been studied many times before, most notably by Smith [17,18], Bogdanova & Ryzhov [19], Pedley & Stephanoff [20], Pedley [16], Pihler-Puzović & Pedley [21]. We follow Pedley & Stephanoff [20] in assuming that the height of the indentation is much greater than the thickness of the boundary layers which form at the walls, but is still small compared with 1, and *σ*=*O*(1). This assumption allows us to proceed to *O*(*ε*^{2}) in the *x*-momentum Navier–Stokes equation to obtain the following equation in the inviscid core:
*y*=0 and *y*=1. It follows that
*P* is eliminated from equations (2.5) and (2.6), the equation for the streamline displacement becomes
*F*(*x*,*t*). They predicted and confirmed experimentally the formation of a wavetrain downstream from the indentation, when the wall was forced with a given frequency. Downstream and upstream from the flexible wall (2.7) reduces to the linearized Korteweg–de Vries (KdV) equation, which supports downstream propagating waves whose group velocity is three times larger than the phase velocity. A common feature of the flow in two-dimensional flexible channels [10,13], these waves were named ‘vorticity waves’, because they owe their existence to the non-zero vorticity gradient in the main flow. For *σ*=*O*(1), the scalings of the model correspond to the lower branch Tollmien–Schlichting (T-S) neutral stability curve, suggesting further that the vorticity waves are a manifestation of the T-S instability. Until recently, it was not established, in general, whether the vorticity waves were a cause or a consequence of the self-excited oscillations in collapsible channels [12,22], because in cases where they had been looked for, vorticity waves occurred together with oscillations (e.g. [10]). However, a recent paper by Xu *et al.* [23] finds an example with oscillations but no vorticity waves, implying that the latter are not an essential cause of the oscillations.

In derivation of equation (2.7), it was implicitly assumed that the boundary layers at the channel walls stay attached at all times and there is no breakaway separation appearing in the flow. However, for some deformations of the wall a small region of reversed flow appears at the upstream corner of the hump as well as at its downstream end, and at a finite time *t*_{cr} the solution to the classical unsteady boundary-layer problem with an adverse pressure gradient breaks down due to the van Dommelen singularity [24]. As soon as breakaway separation of the boundary layer appears (or, in other words, for *t*>*t*_{cr}), (2.7) becomes invalid. Pedley & Stephanoff [20] concluded that if *F*_{x} is *O*(1) and *σ* is *O*(1), then so is *t*_{cr}=*O*((λ^{2}*ε*)^{−1}). As λ^{2}*ε* becomes smaller, the inviscid theory remains valid for a longer time. These estimates were also confirmed by the experimental findings, albeit for the case of prescribed wall oscillations.

Unlike Pedley & Stephanoff [20], we will be solving (2.7) for the elastic wall freely interacting with the fluid. Assuming that the wall is considerably prestretched, the membrane equation for the flexible wall movement is
*T* is the initial longitudinal tension, which is constant up to the required order of accuracy, *P*_{ext} is the pressure external to the channel and *M* represents the mass of the flexible wall. These quantities are related to the dimensional tension *A*,*P*,*F* are functions of time *t* and longitudinal coordinate *x* only.

Pihler-Puzović & Pedley [21] and Kudenatti *et al.* [25] showed that if the pressure is fixed upstream from the membrane region, as well as the flow rate, then problem (2.6)–(2.8) is ill-posed. Therefore, we study these equations only when *P* is fixed downstream. In its essence, the model is inviscid, so if waves propagate in the rigid channels, boundary conditions will be needed to suppress their reflection from the ends. In the experiments, the waves would be attenuated by the presence of viscosity and would (subjected to the rigid channels being sufficiently long) decay before the end of the domain. For the theoretical model formulated here, this implies that the correct boundary conditions need to be given at infinity: far away from the membrane, both upstream and downstream, the flow should be unperturbed. However, as any numerical solution is possible only on a finite domain, we instead apply absorbing boundary conditions (ABCs) at the upstream end of the domain *x*=*b*_{1} (for derivation of the boundary conditions, see appendix A)
*Γ*(*) is the Gamma function. Similarly, if *x*=*b*_{2} is the downstream end of the system, then the ABC for the downstream end has the form

Equations (2.6)–(2.8) subjected to the derived boundary conditions were solved numerically using a combination of Chebyshev collocation points for the spatial discretization with the Crank–Nicolson scheme for the time evolution of the system. The details of the numerical method are presented in appendix B.

## 3. Results

### (a) Solution of the steady problem

In the following, we briefly study the stationary version of the temporal problem (2.6)–(2.8). Although the real interest lies in unsteady results, if the system is stable, in the long-time limit it should converge to a particular solution of the steady problem. The steady equations are invariant with respect to the map 0→1 and 1→0, so we expect symmetry of solutions about *x*=0.5.

#### (i) Large tension, large pressure and small cross-stream pressure gradient limit

The steady problem for *A*=*A*^{s}(*x*), *F*=*F*^{s}(*x*) and *P*=*P*^{s}(*x*) can be solved easily in the asymptotic limit of large tension *T*≫1, large external pressure *P*_{ext}=*TP*_{e} and small cross-stream pressure gradient *σ*≪1. The solution is (see also figure 2)

This solution does not satisfy the steady version of the boundary conditions for the function *A* at the membrane ends (*x*=0 and *x*=1; [26]), because in the limit of *σ*≪1 the highest derivative in the steady version of equation (2.7) is lost and we have to deal with a singular perturbation problem. Therefore, in order to solve the problem completely, we ought to introduce a magnified inner coordinate *X*=*x*/*γ* within an *O*(*γ*) vicinity of each of the membrane ends and consider the appropriate transformed equation. The inner problem does not have an analytical solution, but a relatively simple numerical solution is in good agreement with the steady-state solver (figure 3*a*).

#### (ii) Steady simulations

Other studies of collapsible channels have previously reported a saddle-node bifurcation of the steady channel configurations if the membrane tension *T* was varied, but all other parameters, such as the external pressure *P*_{ext} and the cross-stream pressure gradient parameter *σ*, were kept constant (see [27], where the full interactive boundary layer problem was solved). We find a similar bifurcation structure of the steady solution (figure 4) with one major qualitative difference: in our model the membrane takes only mode-one configurations. Moreover, for *P*_{ext}=0 the only solution to our problem is an unperturbed Poiseuille flow in a parallel-sided channel (or *A*^{s}=*F*^{s}=*P*^{s}=0); the confirmation of this statement comes from solving a perturbation problem to the trivial solution.

For each fixed cross-stream pressure gradient and non-zero external pressure, there exists a critical tension below which the steady problem has no solution. On the other hand, for larger tensions there are two possible steady solutions. The lower branch solution in figure 4*b* corresponds to the less collapsed channel as the tension is increased. An infinite tension gives the limit of a rigid wall, so we anticipate (and confirm later in this paper) that if the system is stable, the unsteady simulations converge to this particular branch. On the other hand, the upper branch in figure 4*c* is never stable and it is clearly unphysical.

Of some interest is the bifurcation point. If *σ* is fixed at a small value, the boundary for existence of the steady solutions represents a parabola (figure 3*b*). Although for higher values of *σ* this dependence is strictly speaking more complicated, the numerical results nevertheless suggest a weak dependence of the boundary position on the value of *σ* (see the different lines in figure 3*b*). When *σ*=0, the steady equations suggest that *A*^{s}(*x*)=−1/2*F*^{s}(*x*). Therefore, the nature of the bifurcation point can be inferred from the problem
*F*^{s}(*x*)=*Tf*(*x*), and equation (3.1) was obtained by combining (2.6) and (2.8). From (3.1) and (3.2), it can be shown that there exists a positive constant *C* such that for

### (b) Unsteady problem

Temporal behaviour of the model was studied by perturbing the steady solution (i) by increasing the external pressure from zero to the constant value to mimic an experimental set-up (*a*), and (iii) by starting the simulations from the steady solution corresponding to a point in parameter space close to the parameter values of interest.

#### (i) *M*=0

In the regime of *M*=0, no fluid–structure instabilities are observed and an ‘inviscid mechanism’, similar to the one reported by Kudenatti *et al.* [25] and Pihler-Puzović & Pedley [21], acts to stabilize the system. When *σ*∼0, the same asymptotic expansion as for the steady solution in §3*a*(i) suggests that problem (2.6)–(2.8) is quasi-stationary at all orders away from the membrane ends (*x*=0 and *x*=1). On the other hand, increasing the influence of the cross-stream pressure gradient by making *σ*∼*O*(1) has little impact on the stability of the system: as the external pressure is increased from zero, the system adjusts to the changes in *P*_{ext} and in the long-time limit always tends to a lower branch stationary solution as in figure 7*a*–*c*. Furthermore, for *M*=0 the solution to the fully coupled nonlinear problem has a distinct algebraic decay to the steady solution (figure 7*d*). This follows from the solution to the linearized KdV equation in the rigid segments [26].

#### (ii) *M*≠0

Adding wall inertia to the model has a dramatic effect on the stability of the system—instabilities are now readily observed, independent of the way the system is perturbed.

If the simulations are started by increasing the external pressure to some finite value of interest, the elastic wall always initially assumes a mode-1 shape in accord with the steady solution in figure 4. However, even when the perturbation to the vector fields decays at first, the system starts oscillating later in time. A typical example is shown in figure 5*a*, where we follow the membrane displacement at *x*=0.5. Note that the solution to the problem when *M*≠0 is exactly superimposed on the decay predicted by the model when *M*=0. By inspecting the corresponding instantaneous shapes *F* of the elastic wall in figure 5*b*, we conclude that the early-time mode-1 deformation is replaced by mode-2 deformations when the system starts oscillating. This instability (its growth rate and frequency) is very robust to the manner in which the simulations are started, so the behaviour of the system can be understood by assuming
*a*(ii). A comparison between this problem and full nonlinear problem (2.6)–(2.8) is presented in figure 6*a* when the spike-like perturbation is introduced to the external pressure. The agreement between the two solutions is excellent if the amplitude of oscillations is not too large. Furthermore, we can assume
*A* is determined by the steady solutions *A*^{s}(*x*), *F*^{s}(*x*) and *P*^{s}(*x*), and the matrix *B* is necessarily singular because of the boundary conditions and the pressure field *P*, which has no explicit time-dependence in the equations (2.6)–(2.8). Once again, the solutions to (3.3) are in good agreement with the nonlinear simulations (figure 6*b*), so the eigenvalue analysis can be used to explore the variation in the system's behaviour with the parameters.

The eigenvalue analysis suggests that the system is neutrally stable only when *σ*=0 (see figure 8*d* in particular). Although the model formally breaks down when *σ*=0, we can still consider the limit when the cross-stream pressure gradient is negligible. In the absence of the pressure gradient, the flow feels the presence of the deformed wall only between *x*=0 and *x*=1, suggesting that the fluid load on the membrane is exactly *P*=−9/2*F*^{2}. Thus, when *σ*=0 the flexible wall behaves in accordance with the equation
*σ* must be non-zero and vorticity waves must be generated in the downstream rigid channel.

Clearly, the wall inertia is crucial for the observed instability. The plots of the fastest growing eigenvalue in figure 8*a*,*b* suggest a complicated dependence on *M*, but, interestingly, *σ*≠0 the frequency of the fastest growing eigenvalues is very close but not equal to a normal mode of the elastic wall *n*, indicated in figure 8*b*, matching the number of extrema in the wall deformation.

From figure 8*c*, it follows that the system stability does not change dramatically if *P*_{ext} is small, and that the instability growth rate stays approximately constant if *P*_{ext}=0. This simplifies the linearized problem significantly, because the corresponding steady solution is unique for all values of *T* and *σ*: it represents an unperturbed Poiseuille flow in the parallel-sided channel with *A*^{s}=*F*^{s}=0.

Therefore, the linearized problem becomes
*F*_{1}/∂*t* for one set of parameters, and to a good approximation it is equal to *D*∂*F*_{1}/∂*t* for constant *D*. Therefore, equation (3.6) essentially behaves like the telegraph equation, for which the solution grows as *n*th mode oscillation is

Finally, although well-posed, this model never leads to constant amplitude oscillations at large times. In other words, the nonlinear terms in equations (2.6)–(2.8) are insufficient to saturate the growth rate of the predicted oscillations. Instead, the instabilities grow until the second branch of steady solutions is encountered and the numerics break down, because the new equilibrium state is exponentially unstable to small amplitude perturbations.

## 4. Discussion of results and conclusion

In this paper, we have studied a rationally derived quasi-one-dimensional model for a collapsible channel that uses the ideas from interactive boundary layer theory. In agreement with the previous findings [7,21,25], the model predicts oscillations, but only if the wall mass is accounted for as well. This result is in qualitative agreement with other studies of flutter in collapsible channels, which predict that the flutter is initially superimposed precisely upon the overall temporal behaviour of the massless system, until the relevant eigenmode is excited [13]. However, unlike the previous work, this model suggests that there is no neutral instability in the system when the wall mass and cross-stream pressure gradient is included.

The model covers instabilities in a long, but finite collapsible channel in the limit of asymptotically high Reynolds number, for which the flow is described by the inviscid equations, but must have some shear. Therefore, the principal features of the resulting flow are independent of the large Reynolds number. The fluid flow creates a loading on the flexible wall, which acts as a negative damping and excites the wall oscillations. The oscillating wall in turn represents a source of the downstream propagating vorticity waves.

Despite the evident simplicity of this model, the ideas of the interactive boundary layer theory can only take us so far. For example, ‘sloshing’-like instabilities predicted in other theoretical work would only come in our model at higher asymptotic orders than the one considered here [8,9,21]. Moreover, the scaling appropriate for T-S instabilities used in this paper is insufficient to capture travelling-wave flutter, which is a structural mode destabilized in the absence of membrane inertia [7]. Similarly, the numerical simulations by Luo & Pedley [10] and Luo *et al.* [11] describe regimes in which the boundary layer thickness *ε* and cross-stream pressure gradient parameter *σ* are of the same order, which is in contrast with the assumptions made here. Nevertheless, many of the assumptions made in this model are essentially the same as in the recent work by Xu *et al.* [12], except that their membrane was even longer than ours, i.e. λ=*O*(*Re*), so that they could reduce the spatial dimension of the system by assuming a parabolic velocity profile everywhere in the channel. This led to a rational choice for both the fluid inertia and the viscous terms. Xu *et al.* [12] found an oscillatory instability even in the absence of wall inertia. As the presence of the viscous terms is the essential difference between the two models, we believe that including suitable ad hoc terms in our equations would also lead to self-excited oscillations when the wall mass is neglected. However, there are no rational ways of improving this model, while remaining in the domain of validity of interactive boundary layer theory, so we proceed no further.

## Funding statement

D.P.-P. thanks Trinity College Cambridge for support through an Overseas Research Scholarship while this work was performed in DAMTP.

## Acknowledgements

The authors are grateful to Prof. Arieh Iserles and Dr Stephen J. Cowley for many helpful discussions.

## Appendix A. The derivation of the boundary conditions

Here, we derive the ABCs, used for the numerical solution of equations (2.6)–(2.8).

In the rigid channels away from the membrane, equation (2.7) is exactly the linearized KdV equation of the form
*et al.* [30], we perform the Laplace transformation on (A 1)
*s*)>0 is the argument in the Laplace-transformed space and *A*. The general solution to (A 2) looks like
_{1}(*s*)=*s*^{1/3}, λ_{2}(*s*)=*ωs*^{1/3}, λ_{1}(*s*)=*ω*^{2} *s*^{1/3}, *c*_{2}(*s*) and *c*_{3}(*s*) must be equal to zero. If *x*=*b*_{1} corresponds to the upstream end of the shortened domain, it follows that
*x*=*b*_{2} is the downstream end of the shortened domain, then a similar procedure leads to the ABC (2.12). A typical validation of the ABCs is shown in figure 10.

It is also useful to consider the leading asymptotic behaviour of (2.9), (2.10) and (2.12) in the long-time limit [26]. This will enable us to compare results of the linear stability analysis with numerical solutions to (2.6)–(2.8), both of which are discussed in §3*b*(ii). Assuming that
*α* is without loss of generality real, then as *α*>0, become
*α*≤0, the boundary conditions are
*b*_{1} and *b*_{2}. If we assume the simplified boundary conditions and consider dependence of the fastest growing eigenvalue on the length of the rigid domains, for sufficiently large values of *b*_{1} and *b*_{2} the growth rate saturates to a constant (figure 11). Throughout §3, the results of the eigenvalue analysis and numerical simulations were presented for *b*_{1}=−10 and *b*_{2}=11.

## Appendix B. Computational method

Strictly speaking, there are three subdomains with different equations and boundary conditions in problem (2.6)–(2.8): the region upstream from the membrane, the membrane region and the region downstream from the membrane. In the upstream and downstream rigid regions, we have to solve (A 1) with (2.9)–(2.10) and (2.12) boundary conditions, respectively. At the junction of the domains, we use conditions of continuity of the streamline displacement function and its derivatives (*A*, *A*_{x} and *A*_{xx}) to match the solution in different domains, so that equivalence between the single and multidomain problems is guaranteed. By taking this approach, we are effectively solving the pressure (2.6) and the membrane (2.8) equations only in the membrane domain (figure 12).

For the spatial discretization of equations (2.6)–(2.8), we used the collocation method based on Chebyshev polynomial expansion (similar to [31]). The equations are therefore satisfied only at the location of nodal points and the grid is denser near the domain ends. Hence, the domain decomposition discussed previously is not only natural, but it also helps to resolve regions around the membrane ends. We map subdomains [*b*_{1},0], [0,1] and [1,*b*_{2}] into the [−1,1] interval, so that the interpolation points represent the extrema of the *N*th-order Chebyshev polynomial and we use Chebyshev collocation differential matrices to approximate *x*-derivatives of functions at the collocation points. The number of collocation points is different in different domains and so are the lengths of the different segments. Therefore, differential matrices are in general different in the three segments.

There are three boundary conditions in each domain, but only two boundary points available. This problem is solved by enforcing one of the boundary conditions implicitly by incorporating it into the differential operators.

For the temporal discretization, the one-step, first-order backward Euler/forward Euler method is used to validate the code when *M*=0. The linear part of the differential operator is evaluated at the new time-level, whereas nonlinearities are computed at the previous time level, so that the linear system that is solved at each time-level has a fixed matrix. This approach speeds up computations significantly. However, as soon as the wall mass is added to the system, the simple implicit Euler is no longer adequate, since it is energy dissipative. Instead, we use the second-order Crank–Nicolson scheme. Consequently, at each time step we apply a Newton–Raphson linearization and solve the linear system for the corrections. All numerical solutions presented in this paper were performed with the time step of Δ*t*=0.001. The problem is solved simultaneously in all three domains for both fluid and structure motions. The resulting matrix is dense since we are using a spectral method for the spatial discretization, but the number of collocation points is relatively small. For *b*_{1}=−10 and *b*_{2}=11,200 collocation points in each of the subdomains *x*∈[*b*_{1},0] and *x*∈[1,*b*_{2}], and 20 collocation points for *x*∈[0,1] was found to be sufficient. Therefore, it is possible to invert the matrices by applying a straightforward LU-decomposition [32]. The corresponding numerical convergence tests were performed by Pihler-Puzović [26].

- Received January 7, 2014.
- Accepted March 20, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.