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.
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 , but are yet to be fully understood in all flow regimes . Moreover, although collapsible channels are inherently global systems with self-excited oscillations often being driven or considerably modified by the presence of finite boundaries , 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  and Stewart et al. . 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 . 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. . 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  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 , 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 , 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 U0, 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: where indicates a dimensional variable, the flow field in the inviscid core, outside the boundary layers formed on the channel walls, is found to be 2.1 2.2 and 2.3where Re1/7≤λ≪Re and Re−2/7≪ε≪1 . Here, u0(y)=6y(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 where the parameter σ=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 , Pedley & Stephanoff , Pedley , Pihler-Puzović & Pedley . We follow Pedley & Stephanoff  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: 2.4and using the no-penetration boundary conditions find the unknown velocity terms from (2.4) in the vicinity of y=0 and y=1. It follows that 2.5and 2.6If P is eliminated from equations (2.5) and (2.6), the equation for the streamline displacement becomes 2.7which was solved previously by Pedley & Stephanoff  for a prescribed 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. ). However, a recent paper by Xu et al.  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 tcr the solution to the classical unsteady boundary-layer problem with an adverse pressure gradient breaks down due to the van Dommelen singularity . As soon as breakaway separation of the boundary layer appears (or, in other words, for t>tcr), (2.7) becomes invalid. Pedley & Stephanoff  concluded that if Fx is O(1) and σ is O(1), then so is tcr=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 , 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 2.8where T is the initial longitudinal tension, which is constant up to the required order of accuracy, Pext is the pressure external to the channel and M represents the mass of the flexible wall. These quantities are related to the dimensional tension , external pressure and membrane mass per unit area via With membrane equation (2.8), we close the fluid–structure interaction problem (2.6) and (2.7). Note that all unknowns A,P,F are functions of time t and longitudinal coordinate x only.
Pihler-Puzović & Pedley  and Kudenatti et al.  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=b1 (for derivation of the boundary conditions, see appendix A) 2.9 and 2.10where 2.11Here Γ(*) is the Gamma function. Similarly, if x=b2 is the downstream end of the system, then the ABC for the downstream end has the form 2.12
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.
(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=As(x), F=Fs(x) and P=Ps(x) can be solved easily in the asymptotic limit of large tension T≫1, large external pressure Pext=TPe 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 ( at x=0 and x=1; ), 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 3a).
(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 Pext and the cross-stream pressure gradient parameter σ, were kept constant (see , 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 Pext=0 the only solution to our problem is an unperturbed Poiseuille flow in a parallel-sided channel (or As=Fs=Ps=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 4b 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 4c 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 3b). 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 3b). When σ=0, the steady equations suggest that As(x)=−1/2Fs(x). Therefore, the nature of the bifurcation point can be inferred from the problem 3.1and 3.2where and Fs(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 there are no solutions to the steady problem [26,28].
(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 ( as in figure 5), (ii) by introducing a short-time spike-like perturbation to the external pressure ( as in figure 6a), and (iii) by starting the simulations from the steady solution corresponding to a point in parameter space close to the parameter values of interest.
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.  and Pihler-Puzović & Pedley , acts to stabilize the system. When σ∼0, the same asymptotic expansion as for the steady solution in §3a(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 Pext and in the long-time limit always tends to a lower branch stationary solution as in figure 7a–c. Furthermore, for M=0 the solution to the fully coupled nonlinear problem has a distinct algebraic decay to the steady solution (figure 7d). This follows from the solution to the linearized KdV equation in the rigid segments .
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 5a, 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 5b, 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 and exploring the solutions to the linearized problem around the steady states computed in §3a(ii). A comparison between this problem and full nonlinear problem (2.6)–(2.8) is presented in figure 6a 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 to obtain a generalized eigenvalue problem 3.3Here, the matrix A is determined by the steady solutions As(x), Fs(x) and Ps(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 6b), 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 8d 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/2F2. Thus, when σ=0 the flexible wall behaves in accordance with the equation 3.4and any perturbation results in the system oscillating around the steady state. Therefore, in order for the instabilities to grow in amplitude, σ 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 8a,b suggest a complicated dependence on M, but, interestingly, is also a singular limit of the model. The reason for this lies in the mechanism of the observed instability. The same mechanism also explains why for σ≠0 the frequency of the fastest growing eigenvalues is very close but not equal to a normal mode of the elastic wall , with the corresponding n, indicated in figure 8b, matching the number of extrema in the wall deformation.
From figure 8c, it follows that the system stability does not change dramatically if Pext is small, and that the instability growth rate stays approximately constant if . In order to infer the nature of the fluid load on the flexible wall and thus explain the mechanism of the instability captured by the model, we consider the limit Pext=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 As=Fs=0.
Therefore, the linearized problem becomes 3.5and 3.6The last term in (3.6) can be interpreted as negative damping. Indeed, in figure 9 we compare the time-evolution of this term against ∂F1/∂t for one set of parameters, and to a good approximation it is equal to D∂F1/∂t for constant D. Therefore, equation (3.6) essentially behaves like the telegraph equation, for which the solution grows as and the frequency of the nth mode oscillation is . Equally important is the last term in (3.5) that excites the fluid, so that there is feedback between the fluid and the elastic wall.
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 . 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 . Similarly, the numerical simulations by Luo & Pedley  and Luo et al.  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. , 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.  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.
D.P.-P. thanks Trinity College Cambridge for support through an Overseas Research Scholarship while this work was performed in DAMTP.
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
In the rigid channels away from the membrane, equation (2.7) is exactly the linearized KdV equation of the form A 1Following Zheng et al. , we perform the Laplace transformation on (A 1) A 2where with ℜ(s)>0 is the argument in the Laplace-transformed space and is the transform of A. The general solution to (A 2) looks like with λ1(s)=s1/3, λ2(s)=ωs1/3, λ1(s)=ω2 s1/3, . Therefore, For the flow field to be unperturbed upstream (at ), c2(s) and c3(s) must be equal to zero. If x=b1 corresponds to the upstream end of the shortened domain, it follows that A 3Finally, the inverse Laplace transformation on (A 3) reveals upstream ABCs (2.9)–(2.10). If x=b2 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 . 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 §3b(ii). Assuming that where α is without loss of generality real, then as the boundary conditions, if α>0, become A 4or, if α≤0, the boundary conditions are A 5These boundary conditions have been validated by solving the generalized eigenvalue problem (3.3) for different values of b1 and b2. 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 b1 and b2 the growth rate saturates to a constant (figure 11). Throughout §3, the results of the eigenvalue analysis and numerical simulations were presented for b1=−10 and b2=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, Ax and Axx) 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 ). 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 [b1,0], [0,1] and [1,b2] into the [−1,1] interval, so that the interpolation points represent the extrema of the Nth-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 b1=−10 and b2=11,200 collocation points in each of the subdomains x∈[b1,0] and x∈[1,b2], 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 . The corresponding numerical convergence tests were performed by Pihler-Puzović .
- Received January 7, 2014.
- Accepted March 20, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.