## Abstract

We present a theoretical description of flow-induced self-excited oscillations in the Starling resistor—a pre-stretched thin-walled elastic tube that is mounted on two rigid tubes and enclosed in a pressure chamber. Assuming that the flow through the elastic tube is driven by imposing the flow rate at the downstream end, we study the development of small-amplitude long-wavelength high-frequency oscillations, combining the results of two previous studies in which we analysed the fluid and solid mechanics of the problem in isolation. We derive a one-dimensional eigenvalue problem for the frequencies and mode shapes of the oscillations, and determine the slow growth or decay of the normal modes by considering the system’s energy budget. We compare the theoretical predictions for the mode shapes, frequencies and growth rates with the results of direct numerical simulations, based on the solution of the three-dimensional Navier–Stokes equations, coupled to the equations of shell theory, and find good agreement between the results. Our results provide the first asymptotic predictions for the onset of self-excited oscillations in three-dimensional collapsible tube flows.

## 1. Introduction

Flow-induced oscillations of fluid-conveying elastic-walled conduits arise in many industrial and biological systems. Examples include pipe flutter, Korotkoff sounds during blood pressure measurement and wheezing during forced exhalation from the pulmonary airways (see the reviews by Heil & Jensen 2003; Grotberg & Jensen 2004). Much current interest in such problems is motivated by oscillations of the vocal cords during phonation (Thomson *et al.* 2005; 2007; Luo *et al.* 2009).

Experimental studies of such systems are usually conducted with a ‘Starling resistor’ set-up (see figure 1), which has become something of a canonical problem in physiological fluid–structure interaction. A length of elastic tubing is connected between two rigid pipes and placed inside a pressure chamber. A flow is driven through the pipes and tubing either by applying a pressure difference between the two ends or by using a volumetric pump at one end. If the transmural pressure (internal minus external) across the tube wall becomes sufficiently negative, the tube will buckle non-axisymmetrically. In this buckled state, small changes in the transmural pressure can induce large changes in the tube shape, resulting in strong fluid–structure interaction. Experiments show that large-amplitude self-excited oscillations tend to develop in the elastic tube when the flow rate exceeds a certain critical value. We refer to Bertram (2003) for a review of relevant experimental studies.

The mechanisms responsible for the development of the self-excited oscillations are still not fully understood. Early theoretical analyses of flow in collapsible tubes were based on spatially one-dimensional models (see the review by Heil & Jensen 2003). To model the complicated three-dimensional wall mechanics, the viscous losses, the effect of flow separation, etc. within the framework of cross-sectionally averaged equations, various ad hoc assumptions are required. Nevertheless, one-dimensional models can still provide useful insights, and capture the qualitative behaviour observed in higher dimensional models, as demonstrated in recent work by Stewart *et al.* (2009). This paper also discusses the link between the widely studied ‘local’ eigenmodes (such as static divergence and travelling wave flutter) of flows in long flexible tubes and channels (Kumaran 1995; Davies & Carpenter 1997; Mandre & Mahadevan 2010) and the ‘global’ eigenmodes (dependent on upstream and downstream boundary conditions) arising in experimental systems such as the Starling resistor.

As a first step towards a more rational theoretical approach to the problem, Jensen & Heil (2003) analysed the development of self-excited oscillations in a two-dimensional channel in which part of one of the rigid walls is replaced by a tensioned membrane. They studied a parameter regime in which the wall stiffness (or pre-stress) is sufficiently large that the wall performs high-frequency oscillations. In this regime, the axial oscillatory fluid flow induced by the wall motion (via conservation of volume) has a two-region structure: an inviscid core region, surrounded by thin Stokes layers near the channel walls. Jensen & Heil (2003) showed that the interaction between the wall-driven oscillatory axial sloshing flow and the mean flow through the tube allows energy to be extracted from the mean flow and provides a mechanism for instability, as described below.

The energy fluxes into and out of the system are the kinetic energy flux through the tube ends, the work done by pressure and viscous forces at the tube ends and at the wall, and the losses due to dissipation. Any net energy flux into the system necessarily alters the total (kinetic and elastic) energy content of the system. If an instability is to grow, the total energy of the system needs to increase, on average, over time, so the net time-averaged energy flux into the system must be positive. By virtue of a non-zero time-mean-square, any oscillatory flow at the upstream end increases the kinetic energy influx, and any oscillatory flow at the downstream end increases the kinetic energy outflux. If the amplitude of the oscillatory flow is greater at the upstream end, then there will be a net increase in kinetic energy flux. If this increase is greater than the additional dissipative (or other) losses, then there is energy available to drive an instability. The mechanism of sloshing instability, identified by Jensen & Heil (2003), therefore relies on an increase in the kinetic energy flux into the system. Other instability mechanisms also exist (Stewart *et al.* 2009), which may instead rely on reductions in dissipative losses, or a decrease in the work done against pressure forces at the tube ends, or possibly a destabilization of internal hydrodynamic modes.

In a theoretical study of the fluid response to prescribed oscillations of the tube wall, Whittaker *et al.* (2010*a*,*b*) confirmed that the essential elements of the sloshing instability mechanism are also present in three-dimensional flow. In a long tube, the additional cross-sectional flows are generally relatively unimportant, and the axial sloshing flows can still dominate the dynamics. Since the sloshing mechanism relies on changes in the cross-sectional area of the pipe, the analysis was restricted to non-circular cross sections. (Small-amplitude deformations of a circular cross section result in area changes that are only quadratic in the displacement amplitude. The induced sloshing flows are then an order of magnitude smaller, resulting in a much smaller influx of kinetic energy and a weaker instability.) Non-circular cross sections are also more representative of the buckled configuration about which oscillations typically develop in experiments.

In this paper, we extend the theoretical work of Whittaker *et al.* (2010*a*) to the case of self-excited oscillations in an elastic tube of initially elliptical cross section. Rather than being imposed externally, the oscillation profile and frequency must now be evaluated as part of the solution. This requires consideration both of the fluid behaviour in response to the wall motion (from Whittaker *et al.* 2010*a*) and of the wall behaviour in response to the pressure in the fluid. For the latter, we make use of the results of Whittaker *et al.* (submitted), in which a ‘tube law’ linking the cross-sectional area of the tube to the transmural pressure is derived from shell theory for an elliptical tube.

In a two-dimensional channel, the restoring mechanisms for wall displacements arise from axial bending and axial tension, and in long-wavelength high-tension limits the bending is typically neglected. For three-dimensional tubes, there are additional bending and stretching effects related to deformations of each cross section. In the long-wavelength limit considered by Whittaker *et al.* (submitted) only azimuthal bending and axial tension are retained at leading order.

This paper is organized as follows. In §2, we define the problem, non-dimensionalize the variables and specify the parameter regime under consideration. In §3, we review the results of Whittaker *et al.* (2010*a*, submitted) for the fluid flow and wall mechanics. In §4, we use these results to derive the leading-order equation governing the normal modes and frequencies of the system. The unsteady three-dimensional fluid–structure interaction reduces to a remarkably simple linear fourth-order ordinary differential equation (4.1). We then use the time-averaged energy budget to determine the growth or decay rates of these modes. We compare our results with direct numerical simulations of the system in §5, and present some concluding remarks in §6.

## 2. Set-up

### (a) Problem description

We consider the set-up depicted in figure 2, with a tube of length *L* and circumference 2*πa*. In its undeformed configuration, the tube has an axially uniform elliptical cross section and is aligned with the coordinate axes (*x**,*y**,*z**), such that the tube axis lies along the *z**-axis, and the major and minor axes of the cross section are aligned with the *x**- and *y**-axes, respectively. Throughout this paper, stars are used to distinguish dimensional quantities from their dimensionless equivalents. The ellipticity is set by a parameter *σ*_{0}, with the major and minor radii being given by and , where
2.1
to obtain the correct circumference . The cross-sectional area in the undeformed configuration is then given by
2.2
where
2.3
is the complete elliptic integral of the second kind.

Sections of the tube adjacent to each end, occupying 0<*z**/*L*<*z*_{1} and *z*_{2}<*z**/*L*<1, are rigid, but for *z*_{1}<*z**/*L*<*z*_{2} the tube wall is made from elastic material which is able to deform in response to the external pressure and the fluid traction. We assume that the tube is subject to a uniform axial pre-stress *F*/(2*πah*) when it is in its undeformed elliptical configuration, owing to an axial tension force *F* applied at the ends of the tube. The wall material is assumed to behave linearly elastically over the range of deformations considered, and has thickness *h*, incremental Young’s modulus *E* and Poisson’s ratio *ν*. From these parameters, we define a bending stiffness
2.4
At *z**=*z*_{1}*L* and *z*_{2}*L*, the elastic section is clamped to the rigid sections, i.e. the position and axial slope of the flexible section are both fixed.

The fluid inside the tube has density *ρ* and viscosity *μ*. We drive the flow by imposing a steady axial volume flux of size at the downstream end. At the upstream end, we fix the pressure *p**=*p**_{up}. (These boundary conditions are chosen to ensure that the instability mechanism described in the introduction is at its most powerful—prescribing the flow rate at the downstream end ensures that no kinetic energy is lost there.) An external pressure *p**_{ext} acts on the outside of the flexible tube.

We consider oscillations of the fluid and elastic wall, with typical amplitude *d*≪*a* and time scale *T*. The key variables used to describe the system will be the axial fluid velocity *w**, the fluid pressure *p**, and the cross-sectional area of the tube *A** at each axial position. For the parameter regime considered here, the transverse velocity components do not appear at leading order.

### (b) Non-dimensionalization

We non-dimensionalize axial lengths on the tube length *L* and transverse lengths on the radial scale *a*, writing
2.5
Times are non-dimensionalized via *t**=*Tt*, where *T* is the expected time scale of the oscillations, derived in equation (2.8) below.

The axial velocity and pressure of the fluid are decomposed into steady and oscillatory components and non-dimensionalized as 2.6 where overbars denote the steady components, and hats the oscillatory perturbations. For the velocity, we non-dimensionalize using the natural scales from the mean flow and wall motion. For the pressure, we use the viscous scale for the steady component, and the inertial scale for the oscillatory component. The steady external pressure is similarly non-dimensionalized as 2.7

The time scale *T* is determined by assuming that the oscillations involve a balance between axial fluid inertia and forces from azimuthal bending of the tube wall. The pressure scale associated with such bending is *Kd*/*a*^{4}. Equating this with the inertial scale *ρL*^{2}*d*/(*aT*^{2}) in equation (2.6)_{2} we obtain
2.8

Finally, energy and energy fluxes are non-dimensionalized on 2.9 based on the kinetic energy and kinetic energy fluxes in the steady flow.

### (c) Dimensionless groups

There are six dimensionless groups in the problem. Three are geometric ratios, corresponding to the wall thickness, tube length and oscillation amplitude,
2.10
There are two additional groups related to the fluid mechanics. Following Whittaker *et al.* (2010*a*) we choose to represent these as a Womersley number *α* and a Strouhal number *St*, defined as
2.11
The Womersley number *α* gives the relative importance of unsteady inertia to viscous effects. Equation (2.11)_{1} shows that, in this problem, *α* is a material parameter. The Strouhal number *St* indicates the relative importance of unsteady to convective inertia. The (steady) Reynolds number is given by
2.12
Finally, there is one additional group related to the tube wall. Following Whittaker *et al.* (submitted), we define the dimensionless axial force
2.13
as the ratio of the restoring forces from tension effects (∼*Fd*/(*aL*^{2})) to those from azimuthal bending (∼*Kd*/*a*^{4}).

We shall consider a parameter regime in which the tube wall is thin, subject to large axial tension, and performs small-amplitude, high-frequency long-wavelength oscillations, so that
2.14
This allows us to use results from Whittaker *et al.* (2010*a*) concerning the fluid flow and from Whittaker *et al.* (submitted) concerning the wall mechanics. (The condition ℓ≫*δ*^{−2/3} was found to be necessary in the latter work, to allow the neglect of forces arising from shearing due to non-axially uniform deformations.)

## 3. Analytical modelling

### (a) Area changes and the steady-state configuration

We assume that the (dimensionless) area *A*(*z*,*t*) of each cross section of the flexible tube varies harmonically in time with dimensionless frequency *ω* and slowly varying dimensionless amplitude Δ(*t*). The combined effects of the steady viscous pressure drop in the fluid and the non-zero external pressure mean that the transmural pressure has a steady, axially varying component. The tube wall responds to this pressure by deforming slightly. The oscillations then take place about this deformed *steady configuration*. We write
3.1
where *A*_{0} is the undeformed cross-sectional area, represents the deformation induced by the steady component of the transmural pressure and describes the (possibly complex) axial mode shape of the oscillatory wall deflection. The scaling for the steady area perturbation is set by ensuring that the ratio *α*^{2}ℓ*St*Δ between the steady and oscillatory perturbations is the same as the ratio between the scalings for the steady and oscillatory pressures in equation (2.6)_{2}. Since *α*^{2}ℓ*St*≫1 and Δ≪1 by assumption, we expect both the oscillatory and steady perturbations to be small. In the rigid sections of the tube (0<*z*<*z*_{1} and *z*_{2}<*z*<1) we impose .

### (b) The fluid response to a moving wall

Whittaker *et al.* (2010*a*) derive an asymptotic description of the fluid flow in response to prescribed periodic oscillations of the tube wall, valid in a high-frequency, long-wavelength, small-amplitude regime, where
3.2
We assume that the growth rate of the free oscillations will be sufficiently slow that the leading-order solutions obtained in Whittaker *et al.* (2010*a*) still hold. (This is justified provided that the additional time derivatives resulting from the time dependence in Δ(*t*) are small compared with those that arise from *e*^{iωt}. This is the case for slowly growing or decaying oscillations that develop when the flow rate is close to the critical value for the onset of self-excited oscillations. This multiple-scales approach will be verified in detail elsewhere.)

Whittaker *et al.* (2010*a*) showed that, at leading order (in terms of the parameters listed in equation (3.2)), the flow comprises steady and oscillatory components, as illustrated in figure 3. The steady flow component is unaffected by the *O*((*α*^{2}ℓ*St*)^{−1},Δ) displacements of the wall described by equation (3.1). It is therefore given by three-dimensional Poiseuille flow in the elliptical undeformed configuration
3.3
The steady pressure is then uniform in each cross section, with a viscous gradient
3.4
Integrating equation (3.4) and applying the boundary condition at *z*=0, we obtain
3.5
Whittaker *et al.* (2010*a*) show that there is also an additional contribution to the steady pressure, induced by nonlinear Reynolds stresses from the oscillatory flow. The contribution to is of *O*(*α*^{2}ℓ*St*Δ^{2}), which may be larger or smaller than the viscous variation (3.5). We omit the details for brevity because this contribution is much smaller for the numerical simulations we describe in §5.

The oscillatory flow component has the same frequency as the wall motion. Because of the high frequency (*α*^{2}≫1), the oscillatory axial velocity has a plug flow profile in the core, with passive Stokes layers adjacent to the tube wall. The magnitude of the plug flow is determined by conservation of mass, given the changing volume of the tube. The pressure is determined from an inertial balance: an axial pressure gradient is required to accelerate the oscillating fluid. The leading-order oscillatory pressure is therefore uniform in each cross section.

As well as evaluating the axial plug flow component, Whittaker *et al.* (2010*a*) solved for the transverse flow, the flow in the Stokes layers, and the first-order corrections to the axial plug core flow. However, the details of these additional components are not needed in the present work.

For the oscillatory components (defined in equation (2.6)), we therefore concentrate on the leading-order axial velocity and pressure, both of which are uniform in each cross section. We write 3.6 and 3.7 where and represent the (complex-valued) axial mode shapes of the oscillatory axial velocity and pressure, respectively.

Since the cross-sectional area *A* is perturbed only slightly from its value *A*_{0} in the undeformed configuration, the leading-order equations arising from conservation of mass and axial momentum are
3.8
Eliminating , we obtain
3.9
This provides a relationship between and , which we shall use in §4 below.

We now consider the energy budget for the fluid, integrated over the volume of the tube, and averaged over one period of the oscillation. For periodic oscillations, Whittaker *et al.* (2010*a*) found the energy budget to be given by
3.10
where is the mean rate of working by the fluid on the wall, is the mean flux of kinetic energy through the tube ends, is the mean rate of working by pressure forces at the tube ends due to the oscillatory flow and is the mean rate of dissipation by the oscillatory flow. (Energy fluxes due to the mean flow cancel out.) Whittaker *et al.* (2010*a*) showed that
3.11

However, the energy budget (3.10) holds only when the amplitude of the oscillations is fixed. To account for the amplitude varying slowly over time, we must add an additional term involving the period-averaged kinetic energy of the oscillatory flow. Changes to over time require a flux of energy to or from the flow, so the energy budget is
3.12
The factor of ℓ*St* enters as a consequence of the energy and energy flux scales (2.9).

The slowly varying amplitude does not affect the leading-order oscillatory plug flow in the core, and hence does not alter the expressions (3.11) for , and . The mean kinetic energy of the oscillatory flow is dominated by the axial plug flow and is given by 3.13 at leading order.

### (c) The wall response to the pressure in the fluid

Whittaker *et al.* (submitted) derive a rational model of the wall mechanics, valid for a thin-walled, initially elliptical tube, loaded by an axially varying transmural pressure, undergoing long-wavelength small-amplitude deformations, for which
3.14
Their analysis provides a ‘tube law’, an expression relating the cross-sectional area *A* to the transmural pressure *p**−*p**_{ext} in the flexible section of the tube,
3.15
where *k*_{0} and *k*_{2} are *O*(1) constants that depend on the ellipticity of the undeformed tube. Values can be found in Whittaker *et al.* (submitted). This tube law expresses the fact that the force exerted on the tube wall by the transmural pressure is balanced by the elastic restoring forces arising from azimuthal bending (proportional to the change in cross-sectional area *A*−*A*_{0}) and those arising from the axial tension acting on axial curvature (proportional to ).

Whittaker *et al.* (submitted) also provide details of the deformation that the tube cross section undergoes as the cross-sectional area changes. Exaggerated deformations of a cross section with increased and decreased cross-sectional areas are shown in figure 3. The full expressions for the deformation are given in Whittaker *et al.* (submitted). Here we just note that control points on the *x*- and *y*-axes are displaced normal to the tube wall with dimensionless amplitudes
3.16

Since the tube law (3.15) is linear, it can be decomposed into steady and oscillatory components. Using equations (2.5), (2.6)_{2}, (3.1) and (3.7), we obtain
3.17

The steady component (3.17)_{1} allows us to compute the perturbation in the steady configuration for a given external pressure , using the expression (3.5) for . The boundary conditions are . The solution in *z*_{1}<*z*<*z*_{2} is straightforward, and we obtain
3.18
where is the steady pressure gradient given in (3.5)_{2}, and
3.19

The oscillatory component (3.17)_{2} of the tube law provides an equation linking and in the flexible section of the tube, which we shall require in §4 below.

We now consider the energy budget for the tube wall. As for the fluid, we integrate over the whole wall and average over a period of the motion. The only source of energy is the work done by the fluid on the wall, , and the only store is elastic energy . (There is no kinetic energy term since we have neglected wall inertia.) Using the scalings (2.9), we therefore have
3.20
The flux is obtained from the energy budget for the fluid problem. The mean elastic energy is evaluated as follows. The instantaneous elastic energy stored in the wall is given by the work required to reach the instantaneous configuration from the undeformed state. As the tube deformation is linear in the pressure, this may be computed by integrating the product of the pressure and area change in each cross section along the length of the tube. Averaging over a period of the oscillation, and subtracting off the energy in the steady configuration, we obtain
3.21
where ^{†} denotes the complex conjugate.

## 4. Solution of the coupled problem

We now combine the results for the fluid and solid mechanics, described in §3, in order to analyse the self-excited oscillations of the coupled system. For this purpose, we combine the leading-order equations governing the oscillatory components of the pressure , axial velocity and cross-sectional area into a single equation. This leading-order description predicts neutrally stable normal modes, which we calculate along with their frequencies. As a more straightforward, but entirely equivalent, alternative to examining the full system at next order, the energy budget is then used to determine the slow growth or decay of the oscillation over time.

### (a) Governing equations and boundary conditions

Eliminating between equations (3.9) and (3.17)_{2}, we obtain a single equation for , valid for *z*_{1}<*z*<*z*_{2},
4.1

In the rigid sections of the tube (0<*z*<*z*_{1} and *z*_{2}<*z*<1), the tube law is replaced by the constraint . From equation (3.9), the governing equation for is then
4.2

When matching between the different regions at *z*=*z*_{1} and *z*_{2}, continuity of pressure and axial volume flux require that and are continuous, i.e.
4.3
Two additional boundary conditions are required on the fourth-order system in the flexible section. Since the flexible wall is clamped to the rigid sections at *z*=*z*_{1} and *z*_{2}, we must have there.^{1} Hence, from equation (3.9), we obtain
4.4
At the two outer ends of the rigid sections, we apply the pressure and volume flux boundary conditions specified in §2. Since the conditions are steady, we require oscillatory perturbations to have zero amplitude there. The flux condition implies since the axial velocity is uniform in the cross section at leading order, and this is converted to a condition on using equation (3.8)_{2}. We therefore obtain
4.5

### (b) The normal modes

We first consider the solutions in the two rigid sections of the tube. Solving equation (4.2) in 0<*z*<*z*_{1} subject to equation (4.5)_{1}, we obtain
4.6
for some constant *G*. Solving equation (4.2) in *z*_{2}<*z*<1 subject to (4.5)_{2} we obtain
4.7
for some other constant *H*.

We now turn to the solution within the flexible section *z*_{1}<*z*<*z*_{2}, where satisfies (4.1). The boundary conditions are (4.4), together with
4.8
which are obtained from applying equation (4.3) to the solutions (4.6) in the rigid sections and eliminating *G* and *H*. The governing equation (4.1) together with the four homogeneous conditions (4.4) and (4.8) define an eigenvalue problem for *ω* in terms of the problem parameters.

The general solution of equation (4.1) is given by 4.9 where 4.10 4.11and 4.12

Substituting equation (4.9) into the boundary conditions (4.4) and (4.8), and eliminating *A*,*B*,*C* and *D*, we obtain the eigenvalue equation
4.13
together with
4.14
Equations (4.13) and (4.14) provide a countable set of solutions for (*m*_{j},*n*_{j}), from which we can determine the normal modes and eigenfrequencies
4.15
Once and *ω* are known, the other variables and can be recovered using equations (3.8)_{2} and (3.9). The solutions for (*m*_{j},*n*_{j}) can only be found numerically.

For illustrative purposes, we consider tubes of ellipticity *σ*_{0}=0.6, with rigid–flexible joints at *z*_{1}=0.1 and *z*_{2}=0.9. For *σ*_{0}=0.6, Whittaker *et al.* (submitted) give *k*_{0}=11.07487 and *k*_{2}=1.70441. Values of *ω*_{j} for such a tube are given in table 1, and some typical normal modes are plotted in figure 4.

Equation (4.15) shows that the non-dimensional frequency depends only on the tube geometry (which also determines the numerical coefficients *k*_{0} and *k*_{2}) and the axial tension parameter, . In dimensional terms, the frequency of the oscillation therefore depends on the tube geometry, the axial tension force *F*, the bending stiffness *K* and the fluid density *ρ* (via equations (2.8) and (2.13)), but not on the mean velocity .

Figure 4 shows the first few normal modes for . As expected, each increase in mode number *j* adds an additional extremum to the axial mode shape. The tension can be varied from 10^{−2} to 10^{2} with little observable variation in the pressure and velocity profiles. The amplitude of the cross-sectional area does show distinct variation in its degree of upstream–downstream symmetry, being more symmetric for higher axial tensions. The tension also has a significant effect on the frequency of the modes, as can be seen in table 1.

### (c) Stability criterion and growth rate

We now determine the slow growth or decay of the normal modes using the energy budget. Eliminating between equations (3.12) and (3.20), we obtain 4.16

Re-writing the expressions (3.13) and (3.21) for and in terms of using equation (3.8) we find that 4.17

The time-averaged energy of the mode is therefore partitioned equally between the kinetic energy in the fluid and elastic energy in the wall. The expressions (3.11) for , and can also be re-written in terms of to yield 4.18

Noting that both of the bracketed expressions in equation (4.16) are proportional to Δ^{2}, we then have
4.19
The amplitude of the oscillation therefore grows or decays exponentially, and we write
4.20
Using equations (4.17) and (4.18), the growth rate *Λ* is
4.21

Examining figure 4, we see that the fundamental mode *p*_{0}(*z*), which has the lowest frequency, also has the largest value of , and this appears to be true in general. Hence, within the parameter regime considered here, the fundamental mode is always the most unstable. It is interesting to note that, even though has little effect on the mode shape, it has a large effect on the frequency and growth rates through the direct factor of in the expression (4.15) for *ω*_{0}.

The critical point at which *Λ*=0 (i.e. neutrally stable oscillations) can be described by a critical Reynolds number *Re*_{c}. From *Re*=*α*^{2}/*St* and equation (4.21), we have
4.22

## 5. Comparison with direct numerical simulations

### (a) The numerical solution

To assess the accuracy of our asymptotic approximation, we performed extensive numerical simulations of the onset of self-excited oscillations in elliptical collapsible tubes. The object-oriented multi-physics finite-element library `oomph-lib` (Heil & Hazel 2007) was used to discretize the three-dimensional unsteady Navier–Stokes equations, coupled to the equations of large-displacement Kirchhoff–Love shell theory. Details of the implementation and validation can be found in Heil & Boyle (2010), where the onset of self-excited oscillations in initially circular cylindrical tubes was studied.

For the computations presented here, we changed the cross sections of the undeformed tube to be elliptical. We prescribed parallel in- and outflow at the upstream and downstream ends of the system, controlled the flow rate at the downstream end, and subjected the fluid at the upstream end to zero axial traction (corresponding to the imposition of zero fluid pressure). The steady external pressure, , was set to a constant value (where is given by equation (3.5)_{2}), equal to the fluid pressure at the midpoint of the collapsible segment with steady Poiseuille flow in the undeformed configuration.

We started the simulations from an initial condition in which the tube wall is in its undeformed configuration and the velocity field is given by steady Poiseuille flow (3.3) for that configuration. To initiate oscillations of a controllable (small) amplitude, we initially increased by a small amount so that the tube wall began to collapse slightly when the simulation was started. We re-set to when the radial displacement of a control point, located on the tube’s minor half axis halfway along the tube (see figure 3*b*), exceeded 0.5 per cent of its initial radius. Following the decay of transients arising from this impulsive start, the wall performed small-amplitude oscillations, of the form depicted in figure 3, about the steady configuration (characterized by ). We extracted the period and the growth or decay rates of the oscillations by fitting the time trace of the control displacement to an exponentially growing or decaying harmonic oscillation.

The thin viscous Stokes layers that develop in the flow fields near the tube wall required much finer spatial resolutions than in the computations of Heil & Boyle (2010). The standard resolution for the computations presented here involved 48 398 degrees of freedom. Selected computations were repeated with an increased spatial resolution involving 75 136 degrees of freedom to assess the mesh independence of the results. Figure 5 shows a representative comparison of the results obtained on the two meshes.

### (b) Numerical results and comparison against theoretical predictions

Figure 5 shows the time trace of the radial displacement of a material control point on the tube wall, for a tube with ellipticity parameter *σ*_{0}=0.6, aspect ratio ℓ=20, dimensionless wall thickness , Poisson’s ratio *ν*=0.49, and dimensionless axial pre-stress . The upstream and downstream rigid tubes each occupy 10 per cent of the tube’s total length, corresponding to *z*_{1}=0.1 and *z*_{2}=0.9. The Womersley number *α* was set to *α*^{2}=50 and results are shown for three Reynolds numbers, *Re*=300, 350 and 400, corresponding to different steady flow rates. Figure 5 shows that, as anticipated, the growth rate increases with *Re*. The oscillations decay for *Re* less than some critical value *Re*_{c} and grow for *Re*>*Re*_{c}, the latter indicating the onset of self-excited oscillations. The period of the oscillations appears to be independent of *Re*, consistent with the asymptotic theory.

The structure of the velocity field is consistent with the assumptions made in our theoretical analysis. Specifically, the oscillatory perturbation to the axial velocity has a large central core region in which the flow resembles plug flow, with thin Stokes layers near the tube walls. Moreover, the pressure distribution, at a given *z* and *t*, is virtually constant across the tube’s cross section, consistent with the predicted flow-field described in §3*b*. Both the velocity profiles and the pressure distribution are strikingly similar to those shown in figs 3 and 4, respectively, in Whittaker *et al.* (2010*b*), in which the wall motion was prescribed.

Figure 6 provides a direct comparison between the numerical and theoretical predictions for the wall shape and flow during the oscillation. All quantities are plotted as functions of the axial position *z*. The numerical results (represented by markers) are shown at one or more of the quarter-period times when their values are close to their maximum amplitude. The theoretical predictions for each field are plotted with lines at the times when their amplitude is maximal. Since the growth over one period is negligible, a single amplitude Δ is chosen to best fit all the numerical data.

Figure 6*a* shows the radial displacement of the tube wall along its major and minor half-axis at an instant when the major and minor half axis are deformed outwards and inwards, respectively. Figure 6*b* shows the axial variation of the tube’s cross-sectional area at two instants when the tube is in a most strongly deformed configuration. The dashed line shows the analytic prediction for the steady state; owing to the steady viscous pressure drop, the tube is slightly expanded in the upstream half and slightly collapsed in the downstream half.

Figure 6*c* illustrates the oscillatory perturbation to the axial velocity on the tube’s centreline at the two instants when the wall moves with maximum velocity. Since a constant flow rate is prescribed at the downstream end, as the wall oscillates, fluid can only be displaced to or from the upstream end of the tube (*z*=0). The magnitude of the axial velocity perturbation therefore increases towards the upstream end. Finally, figure 6*d* compares numerical and theoretical predictions for the pressure. The large markers show the oscillatory perturbation to the fluid pressure at instants when this quantity has its maximum amplitude. The smaller filled markers represent the full pressure (oscillatory component plus steady viscous pressure drop) at intervals a quarter of a period apart. At *t*=1.440 and 1.707 (corresponding to the axial velocity shown in figure 6*c*), the sloshing flow is close to its maximum velocity and its axial acceleration is small. The pressure distribution is therefore dominated by the viscous pressure drop along the tube. However, at *t*=1.305 and 1.571, the wall is in one of its two most strongly deformed configurations and the fluid is subject to its maximum axial acceleration. This acceleration requires large additional axial pressure gradients whose magnitude increases towards the upstream end of the tube where the amplitude of the sloshing flows is greatest. The full pressure distribution reflects the addition of this oscillatory axial pressure gradient to the steady viscous pressure drop.

Overall, the agreement between numerical and theoretical results is excellent, both qualitatively and quantitatively, demonstrating that the flow structure obtained from the numerical solution of the unsteady Navier–Stokes equations and shell theory is indeed captured by our asymptotic analysis.

In figure 7, we compare the frequencies *ω* and growth rates *Λ* for a range of tube lengths with the theoretical predictions (4.15) and (4.21). The theoretical predictions for the frequency of the oscillation are in excellent agreement with the numerical results, with the difference between the two being less than 1.2 per cent for all the cases considered. The theoretical prediction (4.21) consistently over-estimates the growth rate *Λ* of the oscillations, but, qualitatively, the dependence of the growth rate on the Reynolds number and tube length is captured very well.

In figure 8, we plot the critical Reynolds number *Re*_{c} at which the growth rate is zero. This corresponds to neutral oscillations and is therefore the stability boundary for the self-excited oscillations studied here. The points from the numerical simulations are found by interpolating between runs that have close to zero growth rate. The asymptotic prediction is (4.22). Since the asymptotic results over-estimate the growth rates, the predictions for the critical Reynolds numbers are consistently below the numerical results (with the maximum error being about 8 per cent), but qualitatively the dependence on the tube length and the Womersley number is captured very well.

## 6. Discussion

We have developed a theoretical model describing self-excited oscillations in a three-dimensional elastic-walled tube. The model builds on results from separate studies of the fluid and solid mechanics of the system (Whittaker *et al.* 2010*a*, submitted). It is valid for long-wavelength high-frequency small-amplitude oscillations of an initially elliptical thin-walled elastic tube subject to a large axial tension.

The model takes the form of a set of ordinary differential equations and boundary conditions (4.1)–(4.5)_{1} governing the oscillatory component of the pressure as a function of the axial coordinate *z*. Solutions provide predictions for the normal modes (4.9) and associated frequencies (4.15). Consideration of the system’s energy budget predicts that these modes grow or decay exponentially, with rate *Λ* given by equation (4.21). The leading-order oscillations are governed by a balance between axial fluid inertia and the elastic forces in the wall, and decouple from the steady flow. The mode shapes and frequencies are therefore independent of the imposed steady flux, but the flux is crucial in determining the growth rate: the stability criterion (4.22) is expressed in the form of a critical Reynolds number. The fundamental (lowest frequency) mode appears to always have the largest (most positive) growth rate, and so is the most unstable. Viewed as a function of axial position, the fundamental mode has a single maximum in the cross-sectional area variation (equivalently the displacement amplitude). Our predictions for the mode shapes, frequencies and growth rates are all in good agreement with direct numerical simulations.

Our results demonstrate that the instability mechanism first outlined by Jensen & Heil (2003) for two-dimensional channel flow also applies to flow in three-dimensional collapsible tubes. We wish to stress that the instability arises through a genuine interaction between fluid and solid mechanics, rather than through a hydrodynamic instability that is modified by the presence of the elastic tube wall. The fact that our asymptotic predictions agree so well with direct numerical simulations (which would permit the development of oscillations via alternative mechanisms) suggests that, at least in the parameter regime considered here, self-excited oscillations in three-dimensional systems can arise through this mechanism. In other regions of parameter space, self-excited oscillations will almost certainly develop through other mechanisms. For instance, we expect that, at larger Reynolds numbers, the system will become susceptible to additional (fluid–mechanical) instabilities. These either could evolve passively on top of the primary instability studied here or could provide alternative (and, at sufficiently large Reynolds number, possibly more powerful) mechanisms for the onset of self-excited oscillations. We also note that the character of the self-excited oscillations observed in the experiments of Bertram (e.g. Bertram & Tscherry 2006; Bertram *et al.* 2008; Truong & Bertram 2009) is very different from those studied here. In particular, in Bertram’s experiments (where the flow is driven by a prescribed pressure drop) the greatest fluctuations in the flow rate are observed in the downstream rigid tube. This indicates that the oscillations are caused by a different mechanism.

In the instability studied here, the frequency and mode shape of the oscillations are independent of the mean flow rate. This is because the leading-order oscillatory flow is the result of a balance between oscillatory inertia and wall bending, and the oscillatory components decouple from the mean flow. The mean flow does, however, enter into the energy budget and therefore affects the growth rate: the kinetic energy flux and work done at the tube ends are both proportional to the mean flow rate. Higher mean flows provide a greater influx of kinetic energy into the oscillations, while dissipative losses from the oscillatory flow are unaffected.

The critical flow rate at which the oscillations first become unstable is described by a critical Reynolds number, *Re*_{c}, given in equation (4.22). *Re*_{c} is proportional to the tube aspect ratio ℓ=*L*/*a*, so longer tubes are harder to destabilize by the sloshing mechanism. This is because, while the kinetic energy input and dissipation in the Stokes layers are both proportional to the square of the axial velocity (which increases with ℓ), the dissipation has an additional factor of ℓ from the volume of the Stokes layers.

In the present paper, we restricted ourselves to oscillations about an axially uniform elliptical tube, while typical experimental realizations of the Starling resistor tend to employ axisymmetric tubes (which buckle non-axisymmetrically). We note that Heil & Boyle (2010) showed that the instability mechanism analysed here can still operate efficiently in such tubes, though self-excited oscillations tend to develop only from non-axisymmetrically buckled, axially non-uniform equilibrium states. It is unlikely that the onset of self-excited oscillations from such configurations is amenable to a direct analytical description.

The flow boundary conditions employed in this study (flow driven by imposing the flow rate at the downstream end of the system so that no kinetic energy is lost through the outflow) were chosen to ensure that the instability mechanism considered here is at its most powerful, and also for mathematical convenience. This set-up differs from that used in most experimental studies where the flow tends to be driven by an applied pressure drop. Our boundary conditions could be realized experimentally, at least in principle, by driving the flow with a volumetric pump attached to the end of the downstream rigid tube. Conversely, it should be possible to adapt the theory developed here to the case of pressure boundary conditions at both ends of the tube. The resulting analysis would be slightly more complicated, and care would need to be taken to include the effect of the Reynolds-stress-induced adjustment to the mean flow, as described in Whittaker *et al.* (2010*a*,*b*).

The second major constraint for the direct applicability of our theory to existing collapsible tube experiments is that the fluid and wall properties have to fall in the regime given by equation (2.14). In particular, we need to ensure that the wall has negligible inertia relative to the fluid, and that it performs oscillations of sufficiently high frequency. This requires the Womersley number *α* (a material parameter) to be sufficiently large. Using the properties of the rubber tubes employed by Bertram & Elliott (2003) and Truong & Bertram (2009) (*a*=6.5 mm, *h*=1 mm, *ν*=0.5 and *E*=3.15 *MPa*), taking *L*=1.12 *m* as the total length of tubing (rigid plus flexible), and assuming water (*μ*=10^{−3} kg m^{−1} s^{−1} and *ρ*=10^{3} kg m^{−3}) as the working fluid, yields *α*^{2}=42.6. This is comparable to the range of values considered in the present study, suggesting that it should be possible to use existing experimental set-ups to explore the instability mechanism studied here.

## Acknowledgements

The authors would like to acknowledge the financial support of the Engineering and Physical Sciences Research Council to undertake the project of which this work is a part. R.J.W. and S.L.W. are supported by grant EP/D070910/2; and M.H. by grant EP/D670422/1. S.L.W. is also grateful to the EPSRC for funding in the form of an Advanced Research Fellowship.

## Footnotes

↵1 The clamping conditions used in the numerical computations presented in §5 below also impose at

*z*=*z*_{1},*z*_{2}, but our model is not of high enough order to allow this additional boundary condition to be applied. We therefore expect the development of boundary layers near the joints. In Whittaker*et al.*(submitted), it was found that these boundary layers (and also some less obvious shearing effects) have negligible effect on the bulk solution provided ℓ≫*δ*^{−2/3}≫1, and so can be ignored.

- Received December 4, 2009.
- Accepted April 30, 2010.

- © 2010 The Royal Society