TY - JOUR AU - Whittaker, Robert J AB - Summary We examine the effect of wall inertia on the onset of high-frequency self-excited oscillations in flow through an elastic-walled tube. The previous asymptotic model of Whittaker et al. (Proc. Roy. Soc. A466, 2010), for a long-wavelength high-frequency instability in a Starling-resistor set-up, neglected inertia in the tube wall. Here, we extend this model by modifying the ‘tube-law’ for the wall mechanics to include inertial effects. The resulting coupled model for the fluid and solid mechanics is solved to find the normal modes of oscillation for the system, together with their frequencies and growth rates. In the system and parameter regime considered, the addition of wall inertia reduces the oscillation frequency of each mode, however its effect on the stability of the system is not as straightforward. Increasing wall inertia lowers the mean flow rate required for the onset of instability, and is therefore destabilising. However, at higher flow rates the instability growth rate is decreased, and so wall inertia is stabilising here. Overall, the addition of wall inertia decreases the sensitivity of the system to the mean axial flow rate. The theoretical results show good qualitative and reasonable quantitative agreement with direct numerical simulations performed using the oomph-lib framework. 1. Introduction Fluid flow through elastic-walled tubes occurs in many biological systems. In the human body, the cardiovascular, respiratory and digestive systems all use flexible tubes to transport various fluids around the body. As such, the study of flows in elastic tubes is important in understanding the different phenomena that occur in these biological vessels. In the cardiovascular system, the propagation of pulse waves in the arteries is vital for transporting blood to organs and tissues within the body. This is a well-known example, and one-dimensional (1D) models have been developed (1, 2, 3) which are able to adequately explain many properties of the system. The analysis of this problem is helped by the fact that under normal conditions the arteries have a positive transmural (internal minus external) pressure, which allows them to retain a relatively stiff, inflated state. However many blood vessels, such as the veins above the heart and outside the skull have a negative transmural pressure, which causes the vessels to buckle and collapse non-axisymmetrically. These vessels are much more flexible in their buckled state and small changes in fluid pressure can cause large changes in the cross-sectional area. This leads to strong interaction between the fluid and solid mechanics, which induces many interesting phenomena such as flow limitation and self-excited oscillations (4). The collapse of blood vessels and the subsequent effects play an important role in many situations. For example the collapse of blood vessels is believed to be a part of the auto-regulation of blood flow to many internal organs (5, 6), and the external compression of veins in lower limbs is used to prevent deep-vein thrombosis (7, 8). Also, when the brachial artery is compressed by a cuff around the upper arm in blood pressure measurement, flow induced instabilities occur and ‘Korotkoff sounds’ are generated (9, 10). Self-excited oscillations also occur in the airways and are believed to cause a number of respiratory noises. It is thought that flutter instabilities are the cause of respiratory wheezes during forced expiration (11, 12), and controlled flow-induced vibrations of the vocal chords are used in speech production and can be modelled as a collapsible tube system (13). More examples of biological applications can be found in recent review papers (4, 14, 15). Experimental data of fluid flow through elastic-walled tubes is usually obtained using a Starling resistor (16), which is shown in Fig. 1. This comprises an elastic tube which is clamped between two rigid tubes and enclosed in a chamber with a fixed pressure. Fluid is driven through the tubes, either by applying a controlled pressure difference between the ends of the rigid tubes, or by using a volumetric pump to fix a specific flux at one end. If the transmural pressure (internal minus external) across the tube wall becomes sufficiently large and negative, the tube buckles non-axisymmetrically. Once the tube reaches this buckled state, it becomes highly compliant and small changes in the transmural pressure can cause large changes in the tube shape and cross-sectional area. Experiments by Bertram and Tscherry (17) have shown that if the flow exceeds a certain (set-up dependent) critical Reynolds number, the system develops self-excited oscillations. Other relevant experimental studies can be found in the review by Bertram (18). Fig. 1 View largeDownload slide The set-up of the Starling resistor. An elastic tube is clamped between two rigid tubes and is contained in a pressure chamber with fixed pressure $$p^*_{\mathrm{ext}}$$. Flow is driven through the tube using a controlled pressure difference $$p^*_{\mathrm{up}} - p^*_{\mathrm{dn}}$$ between the two ends. Flow can also be driven through the tube by using a volumetric pump to set a particular flux at either end. The pressure $$p^*_{\mathrm{ext}}$$ in the chamber can be modified to control the degree of collapse of the elastic tube Fig. 1 View largeDownload slide The set-up of the Starling resistor. An elastic tube is clamped between two rigid tubes and is contained in a pressure chamber with fixed pressure $$p^*_{\mathrm{ext}}$$. Flow is driven through the tube using a controlled pressure difference $$p^*_{\mathrm{up}} - p^*_{\mathrm{dn}}$$ between the two ends. Flow can also be driven through the tube by using a volumetric pump to set a particular flux at either end. The pressure $$p^*_{\mathrm{ext}}$$ in the chamber can be modified to control the degree of collapse of the elastic tube Early elastic-walled-tube experiments (reviewed in (18)) found a vast array of different types of oscillations spanning a large range of frequencies. However, the mechanisms involved in developing self-excited oscillations are still not fully understood. Many early theoretical analyses of flow in flexible tubes were based on 1D models, which are discussed in the review (4). Typically, these models involve terms describing three quantities: the mass flux or flow rate of the fluid within the tube, the cross-sectional area of the tube and the transmural pressure, all as functions of an axial coordinate and time. Three equations governing the quantities in the models are derived from axial momentum, mass conservation, and a ‘tube law’ for the wall mechanics that relates the transmural pressure to the local cross-sectional area. These kinds of models are still widely used to model networks of collapsible tubes (19, 20). To be able to describe the complicated 3D wall mechanics, viscous dissipation, the effects of flow separation, etc., ad hoc closure assumptions are required. Many of these models are able to capture qualitative effects, such as the onset of self-excited oscillations, that are observed in higher-dimensional models; see for example (21). Pedley (22) introduced a model comprising a 2D channel, where one wall has a segment replaced by a flexible membrane under longitudinal tension. Fluid is driven through the channel by a pressure drop between the two ends of the channel, and the transmural pressure over the membrane determines the initial shape of the membrane. Many 2D models of flow through flexible tubes are based on this system, including (23) and (24). Further examples can be found in the reviews (4) and (15). Using the system constructed by Pedley (22), Jensen and Heil (25) studied a parameter regime where the tension in the wall is large, using a combination of asymptotic analysis and numerical simulation, and determined a simple ‘sloshing’ instability mechanism. In this regime, the system performs high-frequency oscillations, which are governed by a dynamic balance between fluid inertia and large elastic restoring forces. The oscillations of the wall periodically displace fluid from the flexible region of the tube into the rigid regions, which results in axial sloshing flows in the rigid parts of the tube. If the amplitude of these sloshing flows is greater in the upstream rigid section than in its downstream counterpart, then there is a net influx of kinetic energy into the system. If this influx exceeds additional losses, such as viscous dissipation (most of which is found in the boundary layers near the tube walls) and work done by the pressure at the tube ends, then the system can extract energy from the flow to drive an instability. Jensen and Heil (25) used asymptotic techniques to obtain predictions for the frequency and growth rates of modes arising from this instability. They also found their predictions for the critical Reynolds number at which oscillations develop to be in good agreement with numerical simulations. Whittaker et al. (26, 27) showed that the essential components of this sloshing instability mechanism are also present in a 3D flow. However, for efficient extraction of energy from the mean flow to occur, it is necessary that the tube performs oscillations about a non-axisymmetric mean state‡ (26, 28). Hence, this instability is most likely to occur when either the tube’s undeformed cross section is not circular or a tube with an initially axisymmetric cross section has buckled non-axisymmetrically. The case where the tube has a non-circular undeformed cross section was investigated theoretically by Whittaker et al. (29), who combined models for the fluid behaviour in response to prescribed wall motion (26), and the wall behaviour in response to the fluid pressure in the form of a ‘tube law’ linking the transmural pressure with the cross-sectional area (30). The model in (29) is valid for long-wavelength, high-frequency, small-amplitude oscillations of a thin-walled, initially elliptical elastic tube under large axial tension. The predictions made by the model for the mode shapes, frequencies and growth rates of the oscillations, as well as the critical Reynolds number at which oscillations arise, were found to be in good agreement with direct numerical simulations. However, some effects such as wall inertia, axial bending, in-plane shear forces and non-linear effects in the tube wall were neglected to simplify the mathematics within the model. In this article, wall inertia is added to the model of Whittaker et al. (29). This is done by reintroducing the wall inertia terms to the force-balance equations governing the mechanics of the tube wall. Using the force-balance equations, a new ‘tube law’ is derived. Combining this with an asymptotic description of the fluid mechanics of the problem, a complete system for the interaction between the tube wall and the fluid is constructed. Solving this system, countably many oscillatory modes for the instabilities are found and their frequencies are determined. The stability criterion and growth rates of the modes of the oscillations are also determined. We find that the ‘sloshing’ instability mechanism (as described above) still operates. The addition of wall inertia lowers the critical axial flow rate at which the instability first appears, and increases the growth rate of the instability for moderate and low flow rates. However, at higher flow rates, the presence of wall inertia decreases the growth rate of the instability. Hence wall inertia can act as either a destabilising or a stabilising effect, depending on the parameter values. This article is organised as follows. A description of the mathematical set-up used by Whittaker et al. (29) is provided in section 2. In section 3, we derive the model for the oscillatory perturbations, by building on the work of Whittaker et al. (26, 27, 30). As in these works, it is found that the leading-order solutions are composed of a series of neutrally stable normal modes, each with a distinct eigenfrequency. In section 4, we consider the time-averaged energy budget of the system and use this to determine expressions for the (asymptotically slow) growth rates of the oscillatory normal modes. The stability boundary for each mode is expressed in the form of a critical mean-flow Reynolds number. In section 5, we quantify the asymptotic effect that wall inertia has on the frequency and mode shapes of the oscillations, and also the stability and growth rate of each mode. In section 6, the asymptotic predictions are compared with direct numerical simulations, and good agreement is found. Finally, conclusions are presented in section 7. 2. Mathematical set-up 2.1 Problem description We adopt the same set-up as used by Whittaker et al. (29) and consider a tube of length $$L$$ and circumference $$2\pi a$$ with an initially elliptical cross-section, as shown in Fig. 2. The tube axis is aligned with the $$z^{*}$$-axis and the ellipticity of the tube is set by a parameter $$\sigma_{0}$$. Using this parameter, the major and minor radii are given by $$ac\cosh(\sigma_{0})$$ and $$ac\sinh(\sigma_{0})$$, where Fig. 2 View largeDownload slide The set-up used in (29). An initially elliptical elastic-walled tube is clamped between two rigid tubes. Fluid flows from left to right, due to a volume flux condition at the downstream end Fig. 2 View largeDownload slide The set-up used in (29). An initially elliptical elastic-walled tube is clamped between two rigid tubes. Fluid flows from left to right, due to a volume flux condition at the downstream end   \begin{equation} c = \frac{\pi {\mathrm{sech}}(\sigma_{0})}{2 \mathrm{\textbf{E}}({\mathrm{sech}}(\sigma_{0}))}, \qquad \mbox{and} \qquad \mathrm{\textbf{E}}(\phi) = \int_{0}^{\frac{\pi}{2}} (1-\phi^{2}\sin^{2}\vartheta)^{\frac{1}{2}} {\mathrm{d}} \vartheta \end{equation} (2.1a,b) is the complete elliptic integral of the second kind. The constant $$c$$ has been set so that the circumference is $$2\pi a$$. The cross-sectional area in the undeformed state is then   \begin{equation} A^{*}_{0} = \pi a^{2} c^{2} \cosh(\sigma_{0}) \sinh(\sigma_{0}) = a^{2} \frac{\pi^{3} \tanh(\sigma_{0})}{4[\mathrm{\textbf{E}}({\mathrm{sech}}(\sigma_{0}))]^{2}}. \label{eq:undefarea} \end{equation} (2.2) The tube is split into three regions: two rigid sections occupying $$0 < z^{*}/L < z_{1}$$ and $$z_{2} < z^{*}/L < 1$$, and an elastic-walled section within $$z_{1} < z^{*}/L < z_{2}$$ which is clamped onto the rigid tubes at $$z^{*}=z_{1}L, z_{2}L$$. The elastic section is mounted with an axial pre-stress so that in its initial elliptical configuration an axial tension force $$F$$ acts at the ends, giving rise to a uniform axial pre-stress of $$F/(2 \pi ad)$$. In its initial elliptical configuration, the elastic wall has thickness $$d$$ and density $${\rho_{\mathrm{w}}}$$, giving a mass per unit area $$m = {\rho_{\mathrm{w}}} d$$. The elastic wall is susceptible to deformations from forces arising from the transmural (internal minus external) pressure. We assume that it behaves linearly elastically over the range of deformations we consider here, with incremental Young’s modulus $$E$$ and Poisson ratio $$\nu$$. Using these parameters, we define the extensional stiffness $$D$$ and the bending stiffness $$K$$ of the tube wall as   \begin{equation} D = \frac{Ed}{1-\nu^{2}}, \qquad K=\frac{Ed^{3}}{12(1-\nu^{2})}. \end{equation} (2.3a,b) Within the tube, an incompressible Newtonian fluid with density $${\rho_{\mathrm{f}}}$$ and viscosity $$\mu$$ is driven along the tube by the imposition of a steady axial volume flux $$A_{0}^{*}\mathcal{U}$$ at the downstream end, $$z^{*}=L$$. At the upstream end, $$z^{*}=0$$, the pressure is fixed at $$p^{*}=p_{\mathrm{up}}^{*}$$. By prescribing the flow rate at the downstream end, we ensure that no energy is lost to the mean flow there, which in turn, along with the fixed upstream pressure, ensures that the ‘sloshing’ instability mechanism (as described in section 1) is at its most potent. Outside the tube, there is a constant external pressure $$p_{\mathrm{ext}}^{*}$$, which acts on the tube wall. As in the model of (29), we will consider oscillations of the fluid and tube wall with typical timescale $$T$$ and amplitude $$b(t^{*}) \ll a$$, where $$t^{*}$$ is dimensional time. The key variables we will use to describe the system are the fluid pressure $$p^{*}$$, the axial velocity $$w^{*}$$ of the fluid, and the cross-sectional area $$A^{*}$$ of the tube. In the parameter regime we shall be considering, it is found that the pressure and axial velocity are almost uniform in the tube cross-section, and that the transverse velocity components do not appear at leading order. By assuming that oscillations involve a balance between forces from the azimuthal bending of the tube wall and axial fluid inertia, we can estimate the appropriate timescale $$T$$. For this purpose, we equate the inertial pressure scale $${\rho_{\mathrm{f}}} L^{2} b/(aT^{2})$$ associated with oscillations of the fluid with the pressure scale $$Kb/a^{4}$$ associated with azimuthal bending of the tube wall, to obtain   \begin{equation} T = \left(\frac{{\rho_{\mathrm{f}}} a^{3} L^{2}}{K}\right)^{\frac{1}{2}}. \label{eq:nondimt} \end{equation} (2.4) 2.2 Dimensionless groups and parameter regime There are various dimensionless groups in the problem. We have the three geometric ratios   \begin{equation} \delta = \frac{d}{a}, \qquad \ell=\frac{L}{a} \qquad \mathrm{and} \qquad \Delta(t^{*}) = \frac{b(t^{*})}{a}, \end{equation} (2.5a-c) which correspond to the wall thickness, tube length and oscillation amplitude, respectively. We also have two groups related to the fluid mechanics; the Womersley number $$\alpha$$ and the Strouhal number $${\mathit{St}}$$, defined by   \begin{equation} \alpha^{2} = \frac{{\rho_{\mathrm{f}}} a^{2}}{\mu T} = \left(\frac{{\rho_{\mathrm{f}}} K}{a\ell^{2} \mu^{2}}\right)^{\frac{1}{2}} \qquad \mathrm{and} \qquad {\mathit{St}} = \frac{a}{\mathcal{U}T} = \left(\frac{K}{{\rho_{\mathrm{f}}} a^{3} \ell^{2} \mathcal{U}^{2}}\right)^{\frac{1}{2}}. \end{equation} (2.6a,b) The Womersley number represents the relative importance of unsteady inertia to viscous effects and the Strouhal number represents the relative importance of unsteady to convective inertia. Using these, we can define the Reynolds number $$Re$$ as   \begin{equation} Re = \frac{{\rho_{\mathrm{f}}} \mathcal{U} a}{\mu} = \frac{\alpha^{2}}{{\mathit{St}}}. \label{eq:Reynolds} \end{equation} (2.7) Finally, we introduce a dimensionless axial tension $$\mathcal{F}$$ and a dimensionless wall mass $$M$$, defined by   \begin{equation} \mathcal{F} = \frac{aF}{2 \pi K \ell^{2}}, \qquad M = \frac{ma^{4}}{KT^{2}} \equiv \frac{m}{{\rho_{\mathrm{f}}} a \ell^{2}}. \end{equation} (2.8a,b) The dimensionless tension $$\mathcal{F}$$ is the ratio of the restoring forces $$Fb/2\pi aL^{2}$$ from axial tension effects to the restoring forces $$Kb/a^{4}$$ from azimuthal bending. The dimensionless mass $$M$$ is the ratio of wall inertia forces $$mb/T^{2}$$ to the azimuthal bending forces $$Kb/a^{4}$$ or equivalently the forces $${\rho_{\mathrm{f}}} a b \ell^{2}/T^{2}$$ due to the fluid inertia. As in (29), we consider a parameter regime where the tube wall is thin, under a large axial tension and generates small-amplitude, high-frequency, long-wavelength oscillations. We therefore have   \begin{equation} \Delta(t) \ll 1, \quad \alpha \gg 1, \quad \ell {\mathit{St}} \gg 1, \quad \ell \gg 1, \quad \delta \ll 1, \quad \mathrm{and} \quad \mathcal{F} = O(1). \label{eq:parareg} \end{equation} (2.9) 2.3 Non-dimensionalisation Times are scaled on the time-scale $$T$$, transverse lengths on the radial scale $$a$$, and axial lengths on the tube length $$L$$. In particular, we introduce the dimensionless variables   \begin{equation} t = \frac{t^{*}}{T} , \qquad z = \frac{z^{*}}{L}, \qquad A_{0} = \frac{A_{0}^{*}}{a^{2}}, \qquad A = \frac{A^{*}}{a^{2}}. \end{equation} (2.10a-d) It is assumed that the dimensionless area $$A(z,t)$$ varies harmonically in time with dimensionless frequency $$\omega$$ and amplitude $$\Delta(t)=b(t)/a$$. This induces an oscillatory perturbation to the axial flow, which we therefore non-dimensionalise as   \begin{equation} w^{*} = \mathcal{U}\bar{w} + \frac{Lb}{aT}{\mathrm{Re}}\left(\tilde{w}(z)e^{{\mathrm{i}}\omega t}\right) + \dots, \label{eq:nondimw} \end{equation} (2.11) where $$\bar{w}$$ is the steady component and $$\tilde{w}$$ is the (possibly complex) $$z$$-dependent amplitude of the oscillatory component. The scale for the steady flow comes from the imposed flux $$A_0^* \mathcal{U}$$ at the downstream end. The scale for the oscillatory perturbation arises from the continuity equation, given the amplitude and geometry of the wall motion. As discussed in section 3.2 below, the high frequency and large aspect ratio of the system results in the leading-order oscillatory axial velocity $$\tilde{w}$$ being uniform in the cross section of the tube. This axial flow must be driven by an axial pressure gradient. Thus a similar form is required for the fluid pressure $$p^*$$, with both steady and oscillatory components. We write   \begin{equation} p^{*} - p_{\mathrm{up}}^{*} = \frac{\mu L \mathcal{U}}{a^{2}}\bar{p} + \frac{{\rho_{\mathrm{f}}} L^{2}b}{aT^{2}}{\mathrm{Re}}\left(\tilde{p}(z)e^{{\mathrm{i}}\omega t}\right) + \dots, \label{eq:nondimp} \end{equation} (2.12) where $$\bar{p}$$ is the steady component and $$\tilde{p}$$ is the amplitude of the oscillatory component. The steady component has been non-dimensionalised using the viscous scale $$\mu L \mathcal{U}/a^{2}$$, while for the oscillatory component we use the inertial scale $${\rho_{\mathrm{f}}} L^{2} b/ aT^{2} \equiv \Delta K/a^{3}$$. The steady external pressure $$p_{\mathrm{ext}}^{*}$$ and the transmural pressure $${p_{\mathrm{tm}}}^{*}=p^{*}-p^{*}_{\mathrm{ext}}$$ are non-dimensionalised as   \begin{equation} p_{\mathrm{ext}}^{*} - p_{\mathrm{up}}^{*} = \frac{\mu L \mathcal{U}}{a^{2}}\bar{p}_{\mathrm{ext}}, \qquad {p_{\mathrm{tm}}}^{*} = \frac{K}{a^{3}} {p_{\mathrm{tm}}}, \end{equation} (2.13a,b) based on the the viscous scale and the azimuthal bending scale respectively. Combining these with the expression (2.12) for the fluid pressure $$p^{*}$$, we find that the non-dimensional transmural pressure $${p_{\mathrm{tm}}}$$ can be written as   \begin{equation} {p_{\mathrm{tm}}} = \frac{1}{\alpha^{2} \ell {\mathit{St}}}(\bar{p} - \bar{p}_{\mathrm{ext}}) + \Delta(t) {\mathrm{Re}}\left(\tilde{p}(z)e^{{\mathrm{i}}\omega t}\right) + \dots. \label{eq:nondimptm} \end{equation} (2.14) The components of the transmural pressure will act to deform the tube wall, causing both steady and oscillatory changes in the cross-sectional area. The respective scales are determined by balancing the azimuthal bending forces with the transmural pressure, and we obtain   \begin{equation} A(z,t) = A_{0} + \frac{1}{\alpha^{2}\ell {\mathit{St}}}\bar{A}(z) + \Delta(t){\mathrm{Re}}\left(\tilde{A}(z)e^{{\mathrm{i}}\omega t}\right) + \dots, \label{eq:areadec} \end{equation} (2.15) where $$A_{0}$$ is the cross-sectional area in the undeformed state, $$\bar{A}(z)$$ is the change in area due to the steady component of the transmural pressure, and $$\tilde{A}(z)$$ is the axial mode shape of the change in area due to the oscillations of the wall. Energy and energy fluxes are non-dimensionalised using   \begin{equation} {\rho_{\mathrm{f}}} \mathcal{U}^{2}a^{2}L \qquad \mathrm{and} \qquad {\rho_{\mathrm{f}}} \mathcal{U}^{3} a^{3}. \end{equation} (2.16a,b) These are based on the kinetic energy and kinetic energy fluxes in the steady flow. 3. Leading-order model for the oscillatory perturbations 3.1 Tube law for the wall mechanics Starting from the Kirchhoff–Love shell equations (31, 32), Whittaker et al. (30) derived leading-order equilibrium equations for the tube wall in the normal, azimuthal and axial directions. Using the assumptions of small-amplitude, long-wavelength perturbations of a thin-walled tube, Whittaker et al. were able to express the equilibria in terms of a dimensionless azimuthal hoop stress $$\tilde{N}(\tau,z,t)$$, an in-plane shear stress $$\tilde{S}(\tau,z)$$, and the dimensionless displacement functions $$\xi(\tau,z,t)$$ in the normal direction, $$\eta(\tau,z,t)$$ in the azimuthal direction, and $$\zeta_a(z,t)$$ and $$\zeta(\tau,z,t)$$ in the axial direction.§ Here $$z$$ is the dimensionless axial coordinate as defined above, and $$\tau$$ is a dimensionless azimuthal coordinate. To these equilibrium equations, we add the appropriate inertia terms, arising from the components of the wall acceleration, and non-dimensionalised by the transmural pressure scale $$K/a^3$$ from (2.13). The resulting equations (which differ only by the addition of the terms involving $$M$$) are   \begin{eqnarray} \bar{B}\tilde{N} + \frac{\mathcal{F}}{h}\frac{\partial^{2}{\xi}}{\partial{z}^{2}} - \frac{1}{h}\frac{\partial{}}{\partial{\tau}}\left(\frac{1}{h} \frac{\partial{}}{\partial{\tau}}\left(\frac{\beta}{h}\right)\right) + {p_{\mathrm{tm}}} + \frac{M}{h}\ddot{\xi} &=& 0, \label{eq:xiforce1}\\ \end{eqnarray} (3.1)  \begin{eqnarray} \frac{1}{h}\frac{\partial{\tilde{N}}}{\partial{\tau}} +h\frac{{\mathrm{d}} \tilde{S}}{{\mathrm{d}} z} + \frac{\mathcal{F}}{h}\frac{\partial^{2}{\eta}}{\partial{z}^{2}} + \frac{\bar{B}}{h}\frac{\partial{}}{\partial{\tau}}\left(\frac{\beta}{h}\right) + \frac{M}{h}\ddot{\eta} &=& 0, \label{eq:etaforce1}\\ \end{eqnarray} (3.2)  \begin{eqnarray} \frac{1}{\ell}\frac{\partial{}}{\partial{z}}\left(\nu\tilde{N} + 12(1-\nu^{2})\left(\frac{1}{\delta^{2}\ell^{2}}\frac{\partial{\zeta}}{\partial{z}} +\frac{d\zeta_{a}}{dz}\right)\right) + \mathcal{F}\ell\frac{\delta^{2}}{12}\frac{\partial{\tilde{N}}}{\partial{z}} \qquad & & \nonumber\\ +\mathcal{F}\ell(2 - \nu)\left(\frac{1}{\ell^{2}}\frac{\partial^{2}{\zeta}}{\partial{z}^{2}} + \delta^{2}\frac{d^{2}\zeta_{a}}{dz^{2}}\right) + M\ell\left(\frac{1}{\ell^{2}}\ddot{\zeta} + \delta^{2}\ddot{\zeta}_{a}\right) &=& 0, \label{eq:zetaforce1} \end{eqnarray} (3.3) where $$h(\tau) = c(\sinh^{2}(\sigma_{0}) + \sin^{2}\tau)^{\frac{1}{2}}$$ is the coordinate scale factor for $$\tau$$, while   \begin{equation} \bar{B} = -\frac{c^{2}\sinh(2\sigma_{0})}{2h^{3}} \qquad \mbox{and} \qquad \beta = -\frac{2}{c^{2}\sinh(2\sigma_{0})}\frac{\partial{}}{\partial{\tau}}\left(1+\frac{\partial^{2}{}}{\partial{\tau}^{2}} \right)\eta \end{equation} (3.4a,b) are the base-state azimuthal curvature and the azimuthal curvature perturbation respectively. Shell theory gives expressions for the stresses $$\tilde{N}$$ and $$\tilde{S}$$ in terms of the displacements $$\xi$$, $$\eta$$, $$\zeta_a$$ and $$\zeta$$. Hence (3.1)–(3.3) govern the perturbation of the tube wall generated by a given transmural pressure $${p_{\mathrm{tm}}}$$. The method of Whittaker et al. (30) is now followed to reduce these three coupled partial differential equations (PDEs) to a single ordinary differential equation (ODE) relating the cross-sectional area perturbation $$A-A_0$$ to $${p_{\mathrm{tm}}}$$. This involves various asymptotic and other approximations based on the regime (2.9). Omitting the details for brevity, we obtain   \begin{equation} {p_{\mathrm{tm}}}(z,t) = \frac{k_{0}}{A_{0}}(A(z,t)-A_{0}) - \frac{k_{2}\mathcal{F}}{A_{0}}\frac{\partial^{2}{A(z,t)}}{\partial{z}^{2}} + \frac{k_{2}M}{A_{0}}\frac{\partial^{2}{A(z,t)}}{\partial{t}^{2}}, \label{eq:tubelaw} \end{equation} (3.5) where $$k_0(\sigma_0)$$ and $$k_2(\sigma_0)$$ are the numerically determined $$O(1)$$ constants found by (30). This is a ‘tube law’, linking cross-sectional area changes to the transmural pressure. It is is identical to that derived by (30) in the absence of wall inertia, except for the term involving $$M$$. We observe that the coefficient of this new term matches the original term in $$\mathcal{F}$$. This is expected, as the terms in $$M$$ in (3.1) and (3.2) precisely mirror those in $$\mathcal{F}$$, while (3.3) decouples at leading order. As the tube law (3.5) is linear, we can decompose it into steady and oscillatory components. Applying (2.14) and (2.15), we obtain   \begin{align} \bar{p} - \bar{p}_{\mathrm{ext}} &= \frac{k_{0}}{A_{0}}\bar{A}(z) - \frac{k_{2}\mathcal{F}}{A_{0}}\frac{{\mathrm{d}}^{2}\bar{A}(z)}{{\mathrm{d}} z^{2}}, \label{eq:steadp}\\ \end{align} (3.6)  \begin{align} \tilde{p}(z) &= \frac{k_{0}}{A_{0}}\tilde{A}(z) - \frac{k_{2}\mathcal{F}}{A_{0}}\frac{{\mathrm{d}}^{2}\tilde{A}(z)}{{\mathrm{d}} z^{2}} - \frac{k_{2}M\omega^{2}}{A_{0}}\tilde{A}(z) \, , \label{eq:oscip} \end{align} (3.7) which are valid in the flexible section of the tube, $$z_1 < z < z_2$$. 3.2 Fluid mechanics We assume that the fluid flow is governed by the incompressible Navier-Stokes equations. Whittaker et al. (26) showed that the fluid flow in the system can be decomposed into steady and oscillatory components, which decouple at leading order, and only the axial components of the velocity are present. The leading-order steady component $$\bar{w}$$ is the solution of three-dimensional Poiseuille flow in the elliptical undeformed configuration, and is unaltered by the addition of wall inertia. The precise solution for $$\bar{w}$$ is not needed for the following analysis, and so is omitted here. Using a long-wavelength approximation and the property that the oscillatory component of the axial velocity has a high frequency, Whittaker et al. (26) showed, at leading order, that the oscillatory axial velocity $$\tilde{w}$$ has a plug flow profile in the core, with passive thin viscous boundary layers (Stokes layers) near the wall. They also showed that the leading-order oscillatory pressure $$\tilde{p}$$ is uniform in each cross-section. Hence, we have $$\tilde{w} = \tilde{w}(z)$$, $$\tilde{p}=\tilde{p}(z)$$ in the core. The leading-order oscillatory components of the continuity and axial momentum equations adopt the standard forms:   \begin{equation} A_{0}\frac{{\mathrm{d}} \tilde{w}}{{\mathrm{d}} z} + {\mathrm{i}}\omega\tilde{A} = 0, \qquad {\mathrm{i}}\omega\tilde{w} = -\frac{{\mathrm{d}} \tilde{p}}{{\mathrm{d}} z}. \end{equation} (3.8a,b) Eliminating $$\tilde{w}$$ between (3.8a) and (3.8b), we obtain a second relationship between $$\tilde{p}$$ and $$\tilde{A}$$:   \begin{align} \tilde{A} = -\frac{A_{0}}{\omega^{2}}\frac{{\mathrm{d}}^{2}\tilde{p}}{{\mathrm{d}} z^{2}}. \label{eq:Aprel} \end{align} (3.9) 3.3 Boundary and matching conditions The boundary conditions at the tube ends $$z=0,1$$ and the matching conditions at the interfaces at $$z=z_1,z_2$$ between the flexible and rigid sections of the tube are as derived by Whittaker et al. (26). We re-write each of these conditions in terms of the oscillatory pressure perturbation $$\tilde{p}(z)$$, as this is the primary variable we shall be using in section 3.4 below. At the upstream end $$z=0$$, we fix the pressure in the fluid as a constant. Thus the oscillatory component $$\tilde{p}$$ of the pressure must vanish there. At the downstream end $$z=1$$, we fix the axial volume flux. Thus the amplitude $$\tilde{w}$$ of the oscillatory axial plug flow must be zero there. Using (3.8b) to express $$\tilde{w}$$ in terms of $$\tilde{p}$$, we obtain the conditions   \begin{equation} \tilde{p} = 0 \quad \mbox{at} \quad z=0\,, \qquad \frac{{\mathrm{d}}\tilde{p}}{{\mathrm{d}} z} = 0 \quad \mbox{at} \quad z=1\,. \end{equation} (3.10a,b) At $$z=z_{1},z_{2}$$, we require continuity of pressure and axial volume flux in the fluid, and continuity of area perturbation in the wall. The axial volume flux is proportional to $$\tilde{w}$$ and hence $${\mathrm{d}} \tilde{p}/{\mathrm{d}} z$$ from (3.8b). The area perturbation is zero in the rigid sections and is proportional to $${\mathrm{d}}^{2}\tilde{p}/{\mathrm{d}} z^{2}$$ from (3.9) in the flexible section. The full set of matching conditions is therefore   \begin{equation} [\tilde{p}]_{-}^{+} = 0 \, , \quad \left[\frac{{\mathrm{d}}\tilde{p}}{{\mathrm{d}} z}\right]_{-}^{+} = 0\,, \quad \frac{{\mathrm{d}}^{2}\tilde{p}}{{\mathrm{d}} z^{2}} = 0 \quad \mbox{at} \quad z=z_{1},z_{2}\,. \end{equation} (3.11a-c) 3.4 Final model Eliminating $$\tilde{A}$$ between (3.7) and (3.9), we obtain the following equation governing the oscillatory pressure perturbation $$\tilde{p}(z)$$ inside the flexible region of the tube:   \begin{equation} k_{2}\mathcal{F}\frac{{\mathrm{d}}^{4}\tilde{p}}{{\mathrm{d}} z^{4}} +( M\omega^{2}k_{2} - k_{0})\frac{{\mathrm{d}}^{2}\tilde{p}}{{\mathrm{d}} z^{2}} - \omega^{2}\tilde{p} = 0, \quad \mbox{for} \quad z_{1} < z < z_{2}. \label{eq:wiode} \end{equation} (3.12) In the rigid sections of the tube we have $$\tilde{A} = 0$$, and so (3.9) implies   \begin{equation} \frac{{\mathrm{d}}^{2}\tilde{p}}{{\mathrm{d}} z^{2}} = 0, \quad \mbox{for} \quad 00$$ and $$h>0$$ are defined by   \begin{eqnarray} g^{2} &=& \frac{(z_{2}-z_{1})^{2}}{2k_{2}\mathcal{F}}\left[ k_{0} - M\omega^{2}k_{2} + \sqrt{(k_{0} - M\omega^{2}k_{2})^{2}+ 4\omega^{2}k_{2}\mathcal{F}}\right], \label{eq:m1}\\ \end{eqnarray} (3.18)  \begin{eqnarray} h^{2} &=& \frac{(z_{2}-z_{1})^{2}}{2k_{2}\mathcal{F}}\left[M\omega^{2}k_{2} - k_{0} + \sqrt{(k_{0} - M\omega^{2}k_{2})^{2}+ 4\omega^{2}k_{2}\mathcal{F}}\right]. \label{eq:n1} \end{eqnarray} (3.19) Applying the boundary conditions (3.16) to the general solution (3.17) and eliminating the constants $$A$$, $$B$$, $$C$$ and $$D$$, we find the same eigenvalue equation as in (29):   \begin{eqnarray} z_{1}\Big[2gh(1-\cosh{g}\cos{h})+(g^{2}-h^{2})\sinh{g}\sin{h}\Big] \qquad \nonumber\\ - (z_{2}-z_{1})\frac{g^{2}+h^{2}}{gh}\Big[g\sinh{g}\cos{h} + h\cosh{g}\sin{h}\Big] &=& 0, \label{eq:mn2} \end{eqnarray} (3.20) but the expressions above for $$g$$ and $$h$$ are different. We continue to follow the solution method of (29) to solve the eigenvalue system (3.18)–(3.20) for $$\omega$$, but a slight modification is needed as the new terms involving $$M$$ also contain $$\omega$$. We proceed by eliminating $$g$$ and $$\omega$$ from (3.20) to obtain a single equation for $$h$$. First, we obtain an expression for $$\omega$$ in terms of $$g$$ and $$h$$ by considering the product of (3.18) and (3.19):   \begin{equation} \omega^{2} = \frac{g^{2}h^{2}k_{2}\mathcal{F}}{(z_{2}-z_{1})^4}. \label{eq:freqmn} \end{equation} (3.21) Then taking the difference between (3.18) and (3.19), and using (3.21) to eliminate $$\omega$$, we obtain   \begin{equation} g = \left[\frac{\frac{k_{0}(z_{2}-z_{1})^{2}}{k_{2}\mathcal{F}} + h^{2}}{1 + \frac{Mh^{2}k_{2}}{(z_{2}-z_{1})^{2}}}\right]^{\frac{1}{2}}. \label{eq:mn1} \end{equation} (3.22) We now use (3.22) to eliminate $$g$$ from (3.20), giving us an equation to be solved for a single real unknown $$h$$. Solving numerically using Maple, we find countably many solutions for $$h$$. The relationship (3.22) is then used to recover $$g$$ and finally (3.21) is used to find the eigenfrequencies $$\omega$$. We denote the $$n$$th eigenfrequency as $$\omega_{n}$$, with $$\omega_{1}$$ being the fundamental mode. We observe that as $$\omega \rightarrow \infty$$, $$g = O(1)$$ and $$h = O(\omega)$$. Hence, for large $$\omega$$ the eigenvalue equation (3.20) is approximately   \begin{equation} -z_{1} h^{2}\sinh(g)\sin(h)-\frac{(z_{2}-z_{1})h^{2}}{g}\cosh(g)\sin(h)=0. \label{eq:mn2approx} \end{equation} (3.23) Hence we expect to find solutions at $$h \simeq n \pi$$, for large integers $$n$$. This approximation is helpful when computing the numerical solutions, and explains why it is more convenient to solve the system for $$h$$, rather than directly for $$\omega$$. 4. Stability criterion and growth rate The leading-order solution found above comprises oscillatory normal modes that are all predicted to be neutrally stable. At higher orders these modes will be expected to slowly grow or decay in time. We now follow Whittaker et al. (29), who showed that it is possible to derive these growth rates by considering the global energy budget of the system. By subtracting off the mean-flow components, this energy budget can be expressed as   \begin{equation} \frac{{\mathrm{d}}}{{\mathrm{d}} t} \left(\tilde{\mathbb{E}}_{s} + \tilde{\mathbb{E}}_{f}\right) = \frac{1}{\ell {\mathit{St}}}\left(\mathcal{K} - \mathcal{S} - \mathcal{D}\right) . \label{eq:enenbudrel} \end{equation} (4.1) Here, $$\tilde{\mathbb{E}}_{s}$$ is the total dimensionless energy due to oscillations in the tube wall and $$\tilde{\mathbb{E}}_{f}$$ is the dimensionless oscillatory kinetic energy in the fluid, both averaged over a period of the oscillations. On the right-hand side, we have three energy fluxes associated with the oscillatory component of the flow: $$\mathcal{K}$$ is the mean flux of kinetic energy through the ends of the tube due to the oscillatory perturbation, $$\mathcal{S}$$ is the mean rate of working by oscillatory pressure forces that arise at the tube ends, and $$\mathcal{D}$$ is the mean rate of dissipation in the oscillatory viscous Stokes layer adjacent to the tube wall. The energies and fluxes have been non-dimensionalised using the scalings (2.16). By substituting expressions for $$\tilde{\mathbb{E}}_{s}$$, $$\tilde{\mathbb{E}}_{f}$$, $$\mathcal{K}$$, $$\mathcal{S}$$ and $$\mathcal{D}$$ — calculated from the leading-order normal-mode solutions in section 3.5 — into (4.1), we are able to derive an equation for the evolution of the amplitude $$\Delta(t)$$. 4.1 Energy fluxes, fluid energy and wall energy Whittaker et al. (26, 29) showed that $$\mathcal{K}$$, $$\mathcal{S}$$, $$\mathcal{D}$$ and $$\tilde{\mathbb{E}}_{f}$$ are given by   \begin{eqnarray} \mathcal{K} &=& \frac{3}{4} \pi \ell^{2} {\mathit{St}}^{2} \Delta^{2} |\tilde{w}(0)|^{2}, \label{eq:ensou1}\\ \end{eqnarray} (4.2)  \begin{eqnarray} \mathcal{S} &=& \frac{1}{4} \pi \ell^{2} {\mathit{St}}^{2} \Delta^{2} |\tilde{w}(0)|^{2},\label{eq:ensou2} \end{eqnarray} (4.3)  \begin{eqnarray} \mathcal{D} &=& \frac{\pi \ell^{3} {\mathit{St}}^{3} \Delta^{2} (2\omega)^\frac{1}{2}}{2 \alpha} \int_{0}^{1}|\tilde{w}(z)|^{2} {\mathrm{d}} z, \label{eq:ensou3}\\ \end{eqnarray} (4.4)  \begin{eqnarray} \tilde{\mathbb{E}}_{f} &=& \frac{\Delta^{2} {\mathit{St}}^{2} A_{0} \ell^{2}}{4} \int_{0}^{1} |\tilde{w}(z)|^{2} {\mathrm{d}} z. \label{eq:Ef1} \end{eqnarray} (4.5) These expressions concern energies and fluxes in the fluid, and are evaluated purely from the leading-order axial oscillatory fluid flow $$\tilde{w}(z)$$. Hence the expressions are unchanged by the addition of wall inertia. However, a new expression for the energy $$\tilde{\mathbb{E}}_{s}$$ in the wall must be derived here, taking into account not only the elastic energy also but the kinetic energy too. This calculation is detailed in Appendix B, where it is shown that   \begin{equation} \tilde{\mathbb{E}}_{s} = \frac{\Delta^{2} {\mathit{St}}^{2} A_{0} \ell^{2}}{4\omega^{2}} \int_{0}^{1} \left|\tilde{p}'(z)\right|^{2} + 2k_{2}M\left|\tilde{p}''(z)\right|^{2} \, {\mathrm{d}} z. \label{eq:tilEsfinal} \end{equation} (4.6) We note that in the absence of wall inertia ($$M=0$$), (4.6) matches the expression for $$\tilde{\mathbb{E}}_{s}$$ derived in (26, 29). 4.2 Growth rate Now we have expressions for the various energies and fluxes, we can use these along with (4.1) to find an expression for the growth rate of each of the normal modes and thence a stability criterion. From the expressions (4.5) for $$\tilde{\mathbb{E}}_{f}$$ and (4.6) for $$\tilde{\mathbb{E}}_{s}$$, and using (3.8b) to express $$\tilde{w}$$ in terms of $$\tilde{p}$$, we see that   \begin{equation} \tilde{\mathbb{E}}_{s} + \tilde{\mathbb{E}}_{f} = \frac{\Delta^{2} {\mathit{St}}^{2} A_{0} \ell^{2}}{2 \omega^{2}} \int_{0}^{1} |\tilde{p}'(z)|^{2} + k_{2}M|\tilde{p}''(z)|^{2} {\mathrm{d}} z. \label{eq:walfluen} \end{equation} (4.7) From the expressions (4.2)–(4.4) for $$\mathcal{K}$$, $$\mathcal{S}$$ and $$\mathcal{D}$$, and using (3.8b) to express $$\tilde{w}$$ in terms of $$\tilde{p}$$, we have   \begin{equation} \mathcal{K} - \mathcal{S} - \mathcal{D} = \frac{\Delta^{2} {\mathit{St}}^{2} \ell^{2} \pi}{2\omega^{2}} \left[ |\tilde{p}'(0)|^{2} - \frac{\ell {\mathit{St}} (2\omega)^{\frac{1}{2}}}{\alpha} \int_{0}^{1} |\tilde{p}'(z)|^{2} {\mathrm{d}} z \right]. \label{eq:ensoup} \end{equation} (4.8) Substituting (4.7)–(4.8) into (4.1) and evaluating the time derivative on the left-hand side, we obtain   \begin{equation} \frac{{\mathrm{d}} \Delta}{{\mathrm{d}} t} = \left[ \frac{\pi}{2A_{0}}\left(\frac{\frac{|\tilde{p}'(0)|^{2}}{\ell {\mathit{St}}} - \frac{(2\omega)^{\frac{1}{2}}}{\alpha}\displaystyle\int_{0}^{1} |\tilde{p}'(z)|^{2} {\mathrm{d}} z}{\displaystyle\int_{0}^{1} |\tilde{p}'(z)|^{2} + k_{2} M |\tilde{p}''(z)|^{2} {\mathrm{d}} z} \right) \right] \Delta . \end{equation} (4.9) Hence, the amplitude of each normal mode $$\tilde{p}(z)$$ grows or decays exponentially and we may write   \begin{equation} \Delta(t) = \Delta_{0} e^{\Lambda t}, \end{equation} (4.10) where $$\Delta_{0}$$ is the initial amplitude of the oscillations, and the growth rate $$\Lambda$$ is given by   \begin{equation} \Lambda = \frac{\pi}{2A_{0}}\left(\frac{\frac{|\tilde{p}'(0)|^{2}}{\ell {\mathit{St}}} - \frac{(2\omega)^{\frac{1}{2}}}{\alpha}\displaystyle\int_{0}^{1} |\tilde{p}'(z)|^{2} {\mathrm{d}} z}{\displaystyle\int_{0}^{1} |\tilde{p}'(z)|^{2} + k_{2} M |\tilde{p}''(z)|^{2} {\mathrm{d}} z} \right). \label{eq:gro1} \end{equation} (4.11) 4.3 Stability criterion When $$\Lambda = 0$$, neutrally stable oscillations are obtained. We define the critical Reynolds number $${\mathit{Re}}_c$$ to be the Reynolds number (defined by $$Re = \alpha^{2} / {\mathit{St}}$$) at which $$\Lambda = 0$$. (If the other parameters are held constant, then $${\mathit{Re}}$$ is a measure of the mean flow rate through the tube.) Using the expression (4.11) for the growth rate, we find   \begin{equation} {\mathit{Re}}_c = \frac{\alpha \ell (2\omega)^{\frac{1}{2}}}{|\tilde{p}'(0)|^{2}} \int_{0}^{1} |\tilde{p}'(z)|^{2} {\mathrm{d}} z. \label{eq:rec} \end{equation} (4.12) Using this expression for $${\mathit{Re}}_c$$, the growth rate $$\Lambda$$ may be written as   \begin{equation} \Lambda = \frac{ \pi (Re - {\mathit{Re}}_c) |\tilde{p}'(0)|^{2}}{2 A_{0} \ell \alpha^{2} \displaystyle\int_{0}^{1} |\tilde{p}'(z)|^{2} + k_{2}M |\tilde{p}''(z)|^{2} {\mathrm{d}} z}. \label{eq:gro2} \end{equation} (4.13) The growth rate $$\Lambda$$ and critical Reynolds number $${\mathit{Re}}_c$$ depend both on the problem parameters and the particular mode being considered. When $${\mathit{Re}} \leq {\mathit{Re}}_c$$ for every mode, the system will be stable. If $${\mathit{Re}} > {\mathit{Re}}_c$$ for any mode, then the system will be unstable to that oscillatory mode of perturbation. Wall inertia affects $$\Lambda$$ and $${\mathit{Re}}_c$$ both explicitly through the parameter $$M$$, and also implicitly through the mode shapes $$\tilde{p}$$ and frequencies $$\omega$$. 5. Asymptotic predictions 5.1 The effect of wall inertia on the frequencies $$\omega$$ The numerical approach outlined in section 3.5 is used to determine the eigenfrequencies $$\omega_{n}$$ of the model for different values of the axial tension parameter $$\mathcal{F}$$ and wall inertia parameter $$M$$. The tube geometry is kept fixed, with $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$ as in (29). (Using this value for $$\sigma_{0}$$, we have $$k_{0} = 11.07487$$, $$k_{2}=1.70441$$ and $$A_{0} = 2.73060$$.) The results are shown in Table 1. Table 1 Asymptotic predictions for the eigenfrequencies $$\omega_n$$ of the normal modes, for different values of the axial tension parameter $$\mathcal{F}$$ and the wall inertia parameter $$M$$, when $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$ $$\mathcal{F}$$  $$M$$  $$\omega_{1}$$  $$\omega_{2}$$  $$\omega_{3}$$  $$\omega_{4}$$  $$\omega_{5}$$  0.01  0  6.108  19.06  33.80  50.86  70.67  0.01  0.001  6.091  18.60  31.61  44.81  57.98  0.01  0.01  5.948  15.55  21.62  25.74  29.14  0.01  0.1  4.912  7.728  8.563  9.274  10.03  0.01  1  2.389  2.639  2.792  2.985  3.210  0.1  0  6.991  25.73  54.31  94.37  146.5  0.1  0.001  6.970  25.08  50.68  82.96  120.0  0.1  0.01  6.795  20.82  34.41  47.54  60.37  0.1  0.1  5.545  10.24  13.64  17.23  20.87  0.1  1  2.638  3.514  4.481  5.565  6.694  1  0  11.80  59.41  144.3  268.1  431.4  1  0.001  11.76  57.88  134.6  235.6  353.3  1  0.01  11.45  47.96  91.30  135.1  177.8  1  0.1  9.259  23.63  36.29  49.08  61.51  1  1  4.360  8.148  11.96  15.86  19.75  $$\mathcal{F}$$  $$M$$  $$\omega_{1}$$  $$\omega_{2}$$  $$\omega_{3}$$  $$\omega_{4}$$  $$\omega_{5}$$  0.01  0  6.108  19.06  33.80  50.86  70.67  0.01  0.001  6.091  18.60  31.61  44.81  57.98  0.01  0.01  5.948  15.55  21.62  25.74  29.14  0.01  0.1  4.912  7.728  8.563  9.274  10.03  0.01  1  2.389  2.639  2.792  2.985  3.210  0.1  0  6.991  25.73  54.31  94.37  146.5  0.1  0.001  6.970  25.08  50.68  82.96  120.0  0.1  0.01  6.795  20.82  34.41  47.54  60.37  0.1  0.1  5.545  10.24  13.64  17.23  20.87  0.1  1  2.638  3.514  4.481  5.565  6.694  1  0  11.80  59.41  144.3  268.1  431.4  1  0.001  11.76  57.88  134.6  235.6  353.3  1  0.01  11.45  47.96  91.30  135.1  177.8  1  0.1  9.259  23.63  36.29  49.08  61.51  1  1  4.360  8.148  11.96  15.86  19.75  It is seen from the table that the addition of wall inertia significantly reduces the values of the eigenfrequencies $$\omega_{n}$$ for the higher-order modes, even when $$M$$ is small. However for the fundamental modes, the effect of wall inertia is only significant when $$M$$ is larger than around $$0.1$$. Finally, it is noted that for larger values of $$M$$, the eigenfrequencies increase much more slowly as the mode number $$n$$ is increased. 5.2 The effect of wall inertia on $$\tilde{p}$$, $$\tilde{w}$$ and $$\tilde{A}$$ Using the values of the eigenfrequencies $$\omega_{n}$$, we determine the corresponding normal modes of the pressure $$\tilde{p}_{n}$$, the axial velocity $$\tilde{w}_{n}$$ and the area $$\tilde{A}_{n}$$. We obtain $$\tilde{p}_{n}$$ by applying the boundary conditions (3.16) to (3.17)–(3.19). We choose the normalisation   \begin{equation} \int_{0}^{1}|\tilde{p}_{n}'|^{2} {\mathrm{d}} z = 1, \qquad \tilde{p}_{n}'(0)>0 , \end{equation} (5.1a,b) motivated by the form of $${\mathit{Re}}_c$$ in (4.12). Then (3.8b) and (3.9) allow us to find the corresponding $$\tilde{w}_{n}$$ and $$\tilde{A}_{n}$$. In Fig. 3, we show the first five normal modes of $$\tilde{p}_{n}$$, plotted for different values of $$M$$. From the figure, we see little observable change in the fundamental mode for any of the values of $$M$$. However for higher-order modes, there are notable changes as $$M$$ increases. For $$M \gtrsim 0.1$$, rather than being dominated by oscillations in $$z$$ about $$\tilde{p}_{n}=0$$, the eigenfunctions $$\tilde{p}_{n}$$ also have a significant linear component in $$z$$. Fig. 3 View largeDownload slide Asymptotic predictions for the first five normal modes of $$\tilde{p}_{n}$$, for different values of $$M$$, when $$\mathcal{F}=1$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$. The solutions have been normalised such that $$\int_{0}^{1}|\tilde{p}_{n}'|^{2}dz = 1$$ and $$\tilde{p}_{n}'(0) > 0$$ Fig. 3 View largeDownload slide Asymptotic predictions for the first five normal modes of $$\tilde{p}_{n}$$, for different values of $$M$$, when $$\mathcal{F}=1$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$. The solutions have been normalised such that $$\int_{0}^{1}|\tilde{p}_{n}'|^{2}dz = 1$$ and $$\tilde{p}_{n}'(0) > 0$$ Figure 4 shows $${\mathrm{i}}\tilde{w}_{n}$$ for the first five normal modes. Again there is very little difference between the modes for $$M=0,0.01$$. However, in the cases where $$M=0.1$$ and $$M=1$$, the higher-order modes oscillate in $$z$$ about a non-zero value of $${\mathrm{i}}\tilde{w}_{n}$$. Fig. 4 View largeDownload slide Asymptotic predictions for the first five normal modes of $${\mathrm{i}}\tilde{w}_{n}$$, for different values of $$M$$ when $$\mathcal{F}=1$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$. The normalisation is as in Fig. 3. (Note the different vertical scales used on the two lower plots) Fig. 4 View largeDownload slide Asymptotic predictions for the first five normal modes of $${\mathrm{i}}\tilde{w}_{n}$$, for different values of $$M$$ when $$\mathcal{F}=1$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$. The normalisation is as in Fig. 3. (Note the different vertical scales used on the two lower plots) Finally, Fig. 5 shows the first five normal modes of $$\tilde{A}_{n}$$. As before, it is seen that there is not much difference between the modes when $$M=0$$ and $$M=0.01$$. However, in the cases $$M=0.1$$ and $$M=1$$, it can be seen that the fundamental mode $$\tilde{A}_{1}$$ tends towards being symmetric about $$z=0.5$$. It is also noted that unlike the modes for the pressure $$\tilde{p}_{n}$$ and axial velocity $$\tilde{w}_{n}$$, the higher-order modes for the area oscillate in $$z$$ about $$\tilde{A}_{n}=0$$, for all values of $$M$$. Fig. 5 View largeDownload slide Asymptotic predictions for the first five normal modes of $$\tilde{A}_{n}$$, for different values of $$M$$ when $$\mathcal{F}=1$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$. The normalisation is as in Fig. 3. (Note the different vertical scales used on the two lower plots) Fig. 5 View largeDownload slide Asymptotic predictions for the first five normal modes of $$\tilde{A}_{n}$$, for different values of $$M$$ when $$\mathcal{F}=1$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$ and $$\sigma_{0}=0.6$$. The normalisation is as in Fig. 3. (Note the different vertical scales used on the two lower plots) These observations are consistent with the governing equation (3.15) becoming dominated by a balance between the fourth-derivative term and second-derivative term involving $$M$$ as $$M$$ increases. This in turn leads to eigenmodes for $$\tilde{p}$$ that are approximately a linear function plus a sinusoidal oscillation. 5.3 The effect of wall inertia on $${\mathit{Re}}_c$$ and $$\Lambda$$ The critical Reynolds number $${\mathit{Re}}_{c}$$ and growth rate $$\Lambda$$ are found from the normal-mode solutions using (4.12) and (4.13). In Fig. 6, we plot $${\mathit{Re}}_c$$ against $$M$$ for the first four eigenmodes. From the plots, we see that there are significant differences in the behaviour of the modes with odd and even $$n$$ as $$M$$ increases. For the odd modes, $${\mathit{Re}}_c$$ decreases as $$M$$ increases, while for even modes $${\mathit{Re}}_c$$ increases as $$M$$ increases. Thus the odd modes are destabilised by adding wall inertia, while the even modes are stabilised. We also see that, at least for $$0 \le M \le 1$$ the mode with the smallest $${\mathit{Re}}_c$$ and therefore the most unstable, is still the fundamental $$n=1$$ mode. Fig. 6 View largeDownload slide Asymptotic predictions for the critical Reynolds number $${\mathit{Re}}_c$$ as a function of the wall inertia parameter $$M$$ from (4.12), with $$z_{1}=0.1$$, $$z_{2}=0.9$$, and $$\sigma_{0}=0.6$$, for mode numbers $$n\in\{1,2,3,4\}$$. Note the differing behaviour of the odd and even modes: the odd modes are destabilised by increasing wall inertia, while the even modes are stabilised (except at very small $$M$$) Fig. 6 View largeDownload slide Asymptotic predictions for the critical Reynolds number $${\mathit{Re}}_c$$ as a function of the wall inertia parameter $$M$$ from (4.12), with $$z_{1}=0.1$$, $$z_{2}=0.9$$, and $$\sigma_{0}=0.6$$, for mode numbers $$n\in\{1,2,3,4\}$$. Note the differing behaviour of the odd and even modes: the odd modes are destabilised by increasing wall inertia, while the even modes are stabilised (except at very small $$M$$) In Fig. 7, we show plots of $$\partial\Lambda/\partial Re$$ against $$M$$. Since $${\mathit{Re}}$$ can be thought of as a dimensionless measure of the mean flow through the tube, the gradient $$\partial\Lambda/\partial {\mathit{Re}}$$ captures the sensitivity of the growth rate to changes in the mean flow. For all four modes, as $$M$$ increases, the gradient of $$\Lambda$$ for a given $$A_{0}$$, $$\alpha$$ and $$\ell$$ tends to zero. So wall inertia acts as a damping effect, decreasing the sensitivity of the system. However this damping is seen to be much more significant for even modes than for odd modes. Fig. 7 View largeDownload slide Asymptotic predictions for the gradient $$\partial\Lambda/\partial Re$$ of the growth rate as a function of $$M$$, with $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\sigma_{0}=0.6$$, for mode numbers $$n\in\{1,2,3,4\}$$. This demonstrates the sensitivity of the growth rate to the mean-flow Reynolds number. Again note the differing behaviour of the odd and even modes and $$M$$ increases. All the modes become less sensitive to changes in the $${\mathit{Re}}$$ as $$M$$ increases, but the even modes are affected more strongly Fig. 7 View largeDownload slide Asymptotic predictions for the gradient $$\partial\Lambda/\partial Re$$ of the growth rate as a function of $$M$$, with $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\sigma_{0}=0.6$$, for mode numbers $$n\in\{1,2,3,4\}$$. This demonstrates the sensitivity of the growth rate to the mean-flow Reynolds number. Again note the differing behaviour of the odd and even modes and $$M$$ increases. All the modes become less sensitive to changes in the $${\mathit{Re}}$$ as $$M$$ increases, but the even modes are affected more strongly 6. Comparison with direct numerical simulations 6.1 Method for numerical solution To evaluate the accuracy of our asymptotic predictions, we conducted numerical simulations of the onset of self-excited oscillations in elliptical collapsible tubes. In these simulations, the three dimensional unsteady Navier–Stokes equations, coupled to the equations of large-displacement Kirchhoff–Love shell theory, were discretised using the object-oriented multi-physics finite-element library oomph-lib (33). The implementation and validation of these simulations is detailed by Heil and Boyle (34) in their study of the onset of self-excited oscillations in initially circular cylindrical tubes. In the numerical simulations performed here, the cross sections of the undeformed tube were set to be elliptical, and parallel inflow and outflow were imposed at the upstream and downstream ends of the system. The flow rate at the downstream end was controlled and the fluid at the upstream end was subject to zero axial traction (corresponding to setting zero fluid pressure). The steady external pressure $$\bar{p}_{\mathrm{ext}}$$ was chosen to be equal to the fluid pressure that would be obtained at the mid-point of the elastic section under steady Poiseuille flow of the prescribed flux through the tube in its undeformed state. Each simulation was started with an initial condition where the tube wall is in its undeformed configuration with the velocity field within the tube being given by steady Poiseuille flow for this configuration. To initiate small-amplitude oscillations, $$\bar{p}_{\mathrm{ext}}$$ was increased by a small amount so that the tube wall started to collapse slightly at the beginning of the simulation. When the displacement of a control point situated on the tube’s minor half axis halfway along the tube exceeded 0.5% of its initial radius, $$\bar{p}_{\mathrm{ext}}$$ was re-set to its initial value. Once any transients had decayed, the system performed small-amplitude oscillations about the steady configuration, and the period and growth or decay rates of these oscillations were extracted by fitting the time trace of the displacements to an exponentially growing or decaying harmonic oscillation. The standard resolution for the simulations presented here involved 48,398 degrees of freedom. Selected simulations were repeated with an increased spatial resolution with 75,136 degrees of freedom to verify the mesh-independence of the results. 6.2 Comparison of numerical results and asymptotic predictions We compare the asymptotic predictions and the numerical results for   \begin{equation} \sigma_{0} = 0.6, \quad z_{1}=0.1, \quad z_{2} = 0.9, \quad \mathcal{F} = 1, \quad \ell = 15, \quad \alpha^{2} = 50, \label{eq:numericparam} \end{equation} (6.1) where $$\sigma_{0} = 0.6$$ also implies   \begin{equation} A_{0} = 2.73060, \quad k_{0} = 11.07487, \quad k_{2} = 1.70441. \label{eq:numericsigma0param} \end{equation} (6.2) The remaining parameters, $${\mathit{St}}$$ (or equivalently $${\mathit{Re}}$$) and $$M$$ are varied between the different simulations. In Fig. 8, we plot the analytical approximations and numerical calculations for the eigenfrequency $$\omega_{1}$$ of the fundamental mode, first against $$Re$$ for different values of $$M$$ in Fig. 8(a), and then against $$M$$ in Fig. 8(b). It is seen in Fig. 8(a) that as in the asymptotic predictions, the numerical results indicate that varying the Reynolds number $$Re$$ gives negligible change in $$\omega_{1}$$. In both plots, there is excellent agreement between the asymptotic predictions and the numerical simulations. Fig. 8 View largeDownload slide Comparison of the asymptotic predictions and the numerical results for the eigenfrequency $$\omega_{1}$$ of the fundamental mode for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$. In (a), the asymptotic predictions (lines) and numerical results (points) for $$\omega_{1}$$ are plotted against $$Re$$ for $$M=0,0.02,0.2,1,10$$. In (b), the asymptotic predictions (solid line) and the numerical results (points) seen in (a) for $$\omega_{1}$$ are plotted against $$M$$. Note that the numerical results for different values of $$Re$$ and the same $$M$$ are indistinguishable on this scale Fig. 8 View largeDownload slide Comparison of the asymptotic predictions and the numerical results for the eigenfrequency $$\omega_{1}$$ of the fundamental mode for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$. In (a), the asymptotic predictions (lines) and numerical results (points) for $$\omega_{1}$$ are plotted against $$Re$$ for $$M=0,0.02,0.2,1,10$$. In (b), the asymptotic predictions (solid line) and the numerical results (points) seen in (a) for $$\omega_{1}$$ are plotted against $$M$$. Note that the numerical results for different values of $$Re$$ and the same $$M$$ are indistinguishable on this scale In Fig. 9, we compare the asymptotic predictions and numerical results for the variation in $$A-A_{0}$$ along the length of the tube in the cases $$M=0,0.2,1$$. Here, the numerical simulations are plotted using Reynolds numbers close to the analytically predicted values of the critical Reynolds number $${\mathit{Re}}_{c1}$$ of the fundamental mode, so that the growth rate of this mode is small and the amplitude of this mode has little variation throughout the simulation. The values of time $$t$$ were set so that the numerical results for $$A-A_{0}$$ are near their maxima. To obtain the asymptotic prediction for $$A-A_{0}$$, we used the expression (2.15) for $$A$$, together with the asymptotic expression for the steady component $$\bar{A}$$ from (29) and the asymptotic expression for the fundamental oscillatory mode $$\tilde{A}_{1}$$ from section 3.5 above. The amplitude $$\Delta$$ of the oscillatory component was set so that the amplitude of the area variation matches between the asymptotic prediction and numerical results. From Fig. 9, it is seen that as in the asymptotic predictions, the peaks of $$A-A_{0}$$ in the numerical simulations move slightly towards the upstream end of the tube as $$M$$ is increased. We also also observe that there is good agreement between the asymptotic predictions and numerical results for all values of $$M$$. Fig. 9 View largeDownload slide Comparison of theoretical approximations and numerical calculations for the total area variation $$A-A_{0}$$ along the length of the tube, for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$ and $$M\in\{0,0.2,1\}$$. The numerical results (thinner lines) have been calculated for $$Re=248$$, $$t=1.564$$ when $$M=0$$, $$Re=172$$, $$t=3.063$$ when $$M=0.2$$, and $$Re=126$$, $$t=5.449$$ when $$M=1$$. The values of $$Re$$ are chosen to be near the critical Reynolds number $${\mathit{Re}}_{c1}$$ for the fundamental mode so the amplitude of this mode has only slow variation in time, and the values of $$t$$ are chosen to be near the times where the area variation is maximal. The theoretical approximations (thicker lines) are calculated using (2.15), with $$\bar{A}$$ as given in (29) and the $$\tilde{A}_{1}$$ from section 3.5 here. The value of the amplitude $$\Delta$$ of the oscillatory component is set to be $$\Delta=0.114,0.0443,0.0138$$ when $$M=0,0.2,1$$ respectively, to match the amplitude of the area variation between the theoretical approximations and numerical results Fig. 9 View largeDownload slide Comparison of theoretical approximations and numerical calculations for the total area variation $$A-A_{0}$$ along the length of the tube, for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$ and $$M\in\{0,0.2,1\}$$. The numerical results (thinner lines) have been calculated for $$Re=248$$, $$t=1.564$$ when $$M=0$$, $$Re=172$$, $$t=3.063$$ when $$M=0.2$$, and $$Re=126$$, $$t=5.449$$ when $$M=1$$. The values of $$Re$$ are chosen to be near the critical Reynolds number $${\mathit{Re}}_{c1}$$ for the fundamental mode so the amplitude of this mode has only slow variation in time, and the values of $$t$$ are chosen to be near the times where the area variation is maximal. The theoretical approximations (thicker lines) are calculated using (2.15), with $$\bar{A}$$ as given in (29) and the $$\tilde{A}_{1}$$ from section 3.5 here. The value of the amplitude $$\Delta$$ of the oscillatory component is set to be $$\Delta=0.114,0.0443,0.0138$$ when $$M=0,0.2,1$$ respectively, to match the amplitude of the area variation between the theoretical approximations and numerical results The asymptotic prediction (4.12) and numerical calculations for the critical Reynolds number $${\mathit{Re}}_{c1}$$ of the fundamental mode are shown in Fig. 10 as functions of $$M$$. From the figure, we see that the asymptotic prediction captures the qualitative behaviour well but systematically underestimates the value of $${\mathit{Re}}_{c1}$$ by about 13–18%. Fig. 10 View largeDownload slide Comparison of the theoretical approximation (4.12) and numerical calculations for the critical Reynolds number $${\mathit{Re}}_{c1}$$ of the fundamental mode against $$M$$, for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$. The asymptotic prediction is given by the solid line, and the numerical calculations are given by the points Fig. 10 View largeDownload slide Comparison of the theoretical approximation (4.12) and numerical calculations for the critical Reynolds number $${\mathit{Re}}_{c1}$$ of the fundamental mode against $$M$$, for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$. The asymptotic prediction is given by the solid line, and the numerical calculations are given by the points Finally, in Fig. 11, the asymptotic prediction (4.13) and numerical results for the growth rate $$\Lambda_{1}$$ of the fundamental and fastest growing mode are plotted against $$Re$$ for the cases $$M=0,0.2,1,10$$. From the figure, we see that both the theoretical and numerical results vary linearly with the Reynolds number, with the same gradient for each value of $$M$$. Furthermore, although the theoretical results systematically overestimate the value of $$\Lambda_{1}$$, this error decreases with increasing $$M$$ to the point where the theoretical and numerical results are almost indistinguishable in the case of $$M=10$$. Fig. 11 View largeDownload slide Comparison of the theoretical approximation (4.13) and numerical calculations for the growth rate $$\Lambda_{1}$$ of the fundamental mode against $$Re$$, for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$, and $$M\in\{0,0.2,1,10\}$$. The asymptotic predictions are given by the thicker lines, whereas the numerical results are given by the points joined by the thinner lines Fig. 11 View largeDownload slide Comparison of the theoretical approximation (4.13) and numerical calculations for the growth rate $$\Lambda_{1}$$ of the fundamental mode against $$Re$$, for $$\sigma_{0}=0.6$$, $$z_{1}=0.1$$, $$z_{2}=0.9$$, $$\mathcal{F}=1$$, $$\ell = 15$$, $$\alpha^{2}=50$$, and $$M\in\{0,0.2,1,10\}$$. The asymptotic predictions are given by the thicker lines, whereas the numerical results are given by the points joined by the thinner lines Overall, the asymptotic predictions derived in this article are in good agreement with the numerical simulations for the frequency $$\omega_{1}$$ and area variation $$A-A_{0}$$ of the fundamental mode. Although the asymptotic predictions underestimate the value of the critical Reynolds number $${\mathit{Re}}_{c1}$$, and overestimate the value of the growth rate $$\Lambda_{1}$$ of the fundamental mode, they still capture the qualitative behaviour of these quantities. The error in the approximation for $$\Lambda_{1}$$ decreases with increasing wall inertia. Given the number of approximations made to obtain the asymptotic predictions, and the fact that the various ‘small’ parameters are not all that small in the cases considered here, we believe that the discrepancies observed are acceptable, and are confident that the asymptotic system is indeed capturing the essential physics of the full problem. 7. Discussion and conclusions In this article, we have studied the effects of wall inertia on the onset of high-frequency long-wavelength self-excited oscillations in flow through an elastic-walled tube. The previous model of Whittaker et al. (29) has been extended to include inertial resistance to the motion of the elastic wall. The effect of wall inertia in the new model is characterised by a single parameter $$M$$, which is defined in (2.8b) as the ratio of the inertial to azimuthal bending resistance in the wall. We have quantified the effects of $$M$$ on the modes of oscillation and their stability. Asymptotic analysis, based on the limit of high-frequency long-wavelength small-amplitude oscillations in a thin-walled tube, allowed us to reduce the system to a one-dimensional eigenvalue problem (3.15)–(3.16) for the leading-order pressure perturbation amplitude $$\tilde{p}(z)$$. The eigenvalue equation (3.20) was solved numerically to determine the frequencies $$\omega$$, and a countable set of neutrally stable oscillatory normal modes were found. The global energy budget was used to derive the slow growth or decay rates $$\Lambda$$ of each of these modes (a first-order effect, but completely determined by the leading-order solution). As in (29), a critical point was found for each mode to be neutrally stable. This was expressed as a critical mean-flow Reynolds number $${\mathit{Re}}_c$$, with the system being stable for lower Reynolds numbers and unstable for higher Reynolds numbers. The asymptotic expressions for the critical Reynolds number $${\mathit{Re}}_c$$ and the growth rate $$\Lambda$$ were found in (4.12) and (4.13). The fundamental mode with the smallest number of axial oscillations was found to have the lowest-frequency and the highest growth rate. The critical Reynolds number for this mode therefore gave the stability boundary for the system. It is interesting to note that wall inertia does not enter the governing equation (3.15) in the same way that the fluid inertia does. (Compare the $$M\omega^2$$ term associated with wall inertia to the $$\omega^2$$ one associated with fluid inertia.) This effect on the fluid pressure from the wall inertia is directly proportional to the area changes $$\tilde{A}$$ at each $$z$$. However, the effect on the fluid pressure from the fluid inertia is a result of two axial integrals of $$\tilde{A}$$ (one to obtain the axial fluid velocity using (3.8a), and one to integrate the pressure gradient using (3.8b)). Full numerical simulations have also been conducted, for a number of different parameter sets, using the object-oriented multi-physics finite-element library oomph-lib. Our asymptotic predictions compare reasonably well with these numerical simulations when the small parameters in the asymptotics are of size $$O(0.1)$$. The normal modes and frequencies show excellent agreement (Figs. 8 and 9). The growth rates and critical Reynolds numbers agree less well (Figs. 10 and 11), but the qualitative behaviour is captured and the discrepancies are still within acceptable bounds. In both the numerical simulations and the asymptotic predictions, wall inertia (through the parameter $$M$$) can be seen to affect the modes and the stability results in three distinct ways. First it alters the mode shapes $$\tilde{p}(z)$$ through the $$M$$ in the governing equation (3.15). As $$M$$ increases, the complementary function for (3.15) approaches a linear combination of harmonic and linear functions of $$z$$. This increases the symmetry of the modes, as seen in Figs. 3–5. (Physically, the oscillations are governed more by a local balance between wall inertia and axial tension effects, with the symmetry-breaking fluid inertia becoming less important.) As a result of the increased symmetry, the axial sloshing flows increasingly cancel out between adjacent peaks and troughs in the wall displacements $$\tilde{A}$$. For even modes, these cancellations significantly reduce the sloshing flow at the upstream end, and thus weaken the instability. For odd modes, the unmatched peak still gives a sloshing flow at the upstream end, although the magnitude is decreased slightly. The effect enters the expression (4.12) for the critical Reynolds number $${\mathit{Re}}_c$$ through the $$\tilde{p}'(0)$$ factor, where a lower value of $$\tilde{p}'(0)$$ causes a higher $${\mathit{Re}}_c$$, thus giving a more stable mode. As can be seen in Fig. 3, $$\tilde{p}'(0)$$ is significantly reduced for even modes, and hence this is a significant stabilising effect. However, for the fundamental mode (which is odd), this effect only represents a small stabilising effect. Secondly, $$M$$ affects the frequency of the normal modes. As would be expected, the addition of extra inertia reduces the frequency of the normal-mode oscillations. The effect of this in (4.12) is to reduce the critical Reynolds number $${\mathit{Re}}_c$$, thus destabilising the modes. Physically, this effect arises as follows. In (4.2)–(4.4) the kinetic energy flux $$\mathcal{K}$$ and rate of working flux $$\mathcal{S}$$ both contain the same power of $$\omega$$, but the viscous dissipation $$\mathcal{D}$$ has an additional factor of $$\omega^{1/2}$$. (The viscous Stokes layers are $$O(\omega^{1/2})$$ thick, and have $$O(\omega^{-1/2})$$ velocity gradients. The velocity-gradient-squared integrated over the thickness gives the $$\omega^{1/2}$$ factor.) Hence reducing the frequency reduces the dissipative losses relative to the other energy fluxes, thus destabilising the mode. Finally, $$M$$ affects the oscillatory energy contained in each normal mode. This is because the inclusion of inertia in the wall implies the presence of kinetic energy in the wall. This appears as the $$M$$ term in the expression (4.6) for the mean oscillatory energy $$\tilde{\mathbb{E}}_{s}$$ in the wall. (The other term in (4.6) corresponds to the mean elastic potential energy from bending and stretching.) As $$M$$ increases, the energy in the wall for oscillations of a given shape, frequency and amplitude increases. This means that for the amplitude of the mode to increase or decrease by a given amount, more energy must be gained or lost from the system. If the energy fluxes remain fixed then the growth and decay rates of each mode will be damped, so $$|\Lambda|$$ will be smaller. In the expression (4.13) for $$\Lambda$$, this takes effect through the $$M$$ in integral within the denominator. The effect of this can be seen in Fig. 11, where for example at $${\mathit{Re}}=300$$, increasing $$M$$ from $$0$$ to $$10$$ results in a monotonic decrease in the growth rate $$\Lambda_1$$ for both the asymptotic and numerical solutions. Overall, an increase in wall inertia is a destabilising effect, in that it reduces the critical mean-flow Reynolds number for the fundamental mode. This effect is brought about primarily through the lowering of the frequency of the mode, and the resultant decrease in dissipative losses in the fluid.∥ However, the increase in energy in the oscillations caused by the kinetic energy of the wall means that the growth rates are damped by increasing amounts as more wall inertia is added. At higher Reynolds numbers this effect dominates, and the growth rates are reduced by an increase in wall inertia. In this regime wall inertia acts as a stabilising effect. To assess the range of values of $$M$$ in physiological applications and in collapsible tube experiments we rewrite (2.8b) as   \begin{equation} M =\left(\frac{{\rho_{\mathrm{w}}}}{{\rho_{\mathrm{f}}}}\right) \left(\frac{d}{a}\right) \left(\frac{\vphantom{h}a}{L}\right)^2. \end{equation} (7.1) This shows that large values of $$M$$ arise in scenarios where the solid-to-fluid density ratio is large and the tube walls are relatively thick. An increase in tube length decreases $$M$$, though it is important to stress that in physiological applications where collapsible vessels are not mounted on rigid upstream and downstream supports, $$L$$ should be interpreted as the length over which the vessel collapses; similar considerations apply to experiments with a Starling resistor where collapse often occurs only near the downstream end of a relatively long collapsible tube. Furthermore, we reiterate that our theory only applies directly to situations where the volume flux is controlled at the downstream end, and fluctuations in the flow rate (induced by the changes in the volume of the oscillating collapsible segment) occur at the upstream part of the system. Almost all experiments in the literature were driven by an applied pressure drop and typically the fluctuations in the flow rate downstream of the collapsible segment are then found to be larger than those upstream. The lengths and wall thicknesses of collapsible tubes used in experiments with the Starling resistor reported in the literature span a wide range. Those used by Katz et al. in (35) ($$\ell = 7$$, $$d/a \approx 1/4$$) were relatively short while those used by Bertram and Castles in (36) ($$\ell = 17$$, $$d/a \approx 1/3$$) are more typical. Both sets of tubes are rather thick by the standards of thin shell theory. If flow-induced collapse occurs in the pulmonary airways, the collapse is likely to occur between successive bifurcations which are typically separated by a few vessel diameters. $$M$$ can therefore span a wide range of values, from $$M \approx 1 \times 10^{-3}$$ in (36) for experiments with silicon rubber tubes conveying water ($${\rho_{\mathrm{w}}}/{\rho_{\mathrm{f}}} \approx 1$$) to $$M \approx 3.3$$ for a pulmonary airway of thickness $$d/a = 1/10$$ conveying air ($${\rho_{\mathrm{w}}}/{\rho_{\mathrm{f}}} \approx 816$$) and collapsing over a length of $$\ell = 5$$. The work presented here is not directly comparable with previous theoretical studies of the effect of wall inertia on oscillations in collapsible tubes and channels. Previous studies have concentrated on different regions of parameter space (often lower axial tension and/or a significantly collapsed/buckled base state) and hence have explored different instability mechanisms that are affected by wall inertia in different ways. For example, in the case of channel flow with a lower-tension elastic membrane, Luo and Pedley (37) found that the addition of wall inertia induces an additional higher-frequency ‘flutter’ instability superposed on the original lower-frequency oscillation that occurs in the same setup with a mass-less wall. In the present work, the addition of wall inertia solely acts as a perturbation to the original high-frequency oscillatory modes. Acknowledgements M.C.W. gratefully acknowledges support from the University of East Anglia in the form of a funded Ph.D. Studentship. APPENDIX A Proof that the eigenfrequencies are real In this Appendix, we prove that the eigenfrequencies $$\omega$$ of the system (3.10)–(3.13) always take non-zero real values. We start with the governing ODE (3.12) for $$\tilde{p}$$ in the flexible region ($$z_{1} < z < z_{2}$$) of the tube, multiply it by the complex conjugate of $$\tilde{p}''$$ and integrate between $$z_1$$ and $$z_2$$:   \begin{equation} \int_{z_1}^{z_2} \left\{ k_{2}\mathcal{F}\tilde{p}''''\tilde{p}''^{\dagger} +( M\omega^{2}k_{2} - k_{0})|\tilde{p}''|^{2} - \omega^{2}\tilde{p}''^{\dagger}\tilde{p} \right\} {\mathrm{d}} z = 0, \end{equation} (A.1) where $$'$$ denotes a derivative with respect to $$z$$ and $$^{\dagger}$$ represents the complex conjugate. Using integration by parts on the first and last terms in the integrand, we obtain   \begin{equation} - k_{2}\mathcal{F} \int_{z_{1}}^{z_{2}}|\tilde{p}'''|^{2}\, {\mathrm{d}} z + (M\omega^{2}k_{2}-k_{0})\int_{z_{1}}^{z_{2}}|\tilde{p}''|^{2}\, {\mathrm{d}} z - \omega^{2}\left(\Bigl[\tilde{p}'^{\dagger}\tilde{p}\Bigr]_{z_1}^{z_{2}} - \int_{z_{1}}^{z_{2}}|\tilde{p}'|^{2} \, {\mathrm{d}} z\right) = 0, \label{eq:refreqprf1} \end{equation} (A.2) where the boundary contributions from the first term vanish because $$\tilde{p}''=0$$ at $$z=z_{1},z_{2}$$ from (3.11c). We now perform an analogous procedure with the governing equation (3.13) in the rigid parts of the tube. Taking the complex conjugate of (3.13), multiplying by $$\tilde{p}$$, and integrating over each of the rigid sections, we have   \begin{equation} \int_0^{z_1} \tilde{p}''^{\dagger}\tilde{p} \, {\mathrm{d}} z = 0, \qquad \int_{z_2}^1 \tilde{p}''^{\dagger}\tilde{p} \, {\mathrm{d}} z = 0. \end{equation} (A.3a,b) Using integration by parts on these two integrals, we then obtain   \begin{equation} \omega^{2}\int_{0}^{z_{1}}|\tilde{p}'|^{2}\, {\mathrm{d}} z - \omega^{2}\Bigl[\tilde{p}'^{\dagger}\tilde{p}\Bigr]_{0}^{z_{1}} = 0, \qquad \omega^{2}\int_{z_{2}}^{1}|\tilde{p}'|^{2}\, {\mathrm{d}} z - \omega^{2}\Bigl[\tilde{p}'^{\dagger}\tilde{p}\Bigr]_{z_{2}}^{1} = 0. \end{equation} (A.4a,b) We now sum equations (A.2), (A.4a) and (A.4b). The boundary terms at $$z=z_1,z_2$$ cancel by virtue of the continuity conditions (3.11). Noting that $$\tilde{p}'' = \tilde{p}''' = 0$$ for $$z \in (0, z_{1})$$ and $$z \in (z_{2}, 1)$$, we deduce   \begin{equation} (M\omega^{2}k_{2}-k_{0})\int_{0}^{1}|\tilde{p}''|^{2}\, {\mathrm{d}} z - k_{2}\mathcal{F}\int_{0}^{1}|\tilde{p}'''|^{2}\, {\mathrm{d}} z - \omega^{2}\Bigl[\tilde{p}'^{\dagger}\tilde{p}\Bigr]_{0}^{1} + \omega^{2}\int_{0}^{1}|\tilde{p}'|^{2} \, {\mathrm{d}} z = 0. \label{eq:refreqprf3} \end{equation} (A.5) Finally, we know that $$\tilde{p} = 0$$ at $$z=0$$ from (3.10a), and $$\tilde{p}'=0$$ at $$z=1$$ from (3.10a), so the remaining boundary terms vanish. Rearranging (A.5), we can then write   \begin{equation} \omega^{2} = \frac{\displaystyle\int_{0}^{1}k_{2}\mathcal{F}|\tilde{p}'''|^{2} + k_{0}|\tilde{p}''|^{2} \, {\mathrm{d}} z}{\displaystyle\int_{0}^{1}|\tilde{p}'|^{2} + k_{2}M|\tilde{p}''|^{2}\, {\mathrm{d}} z}. \label{eq:omega-quotient} \end{equation} (A.6) The constants $$k_0$$ and $$k_2$$ are both strictly positive, while $$\mathcal{F}$$ and $$M$$ are non-negative. Hence for any non-trivial solution $$\tilde{p}$$, the right-hand side of (A.6) is real and strictly positive. Hence, $$\omega$$ must be real and non-zero. APPENDIX B In this Appendix, we derive an expression for $$\tilde{\mathbb{E}}_{s}$$, the period-averaged dimensionless oscillatory energy in the tube wall, including both the elastic and kinetic contributions. B.1 Rate of working on the tube wall We start from the following (dimensional) expression for the rate of working on the tube wall due to pressure forces   \begin{equation} \frac{{\mathrm{d}} E^{*}_{s}}{{\mathrm{d}} t^{*}} = \iint_{\mathrm{Tube\;Wall}} {p_{\mathrm{tm}}}^{*} \frac{\partial \mathbf{r}}{\partial t^{*}} \cdot \hat{\mathbf{{n}}} \, {\mathrm{d}} S. \label{eq:dimrow} \end{equation} (B.1) Here, $$E_{s}^{*}$$ is the total dimensional energy in the tube wall as a result of work by the transmural pressure $${p_{\mathrm{tm}}}^{*}$$, $${\mathrm{d}} S$$ is an element of the mid-plane of the tube wall, $$t^{*}$$ is once again dimensional time, $$\mathbf{r}$$ is the position of the wall mid-plane in the deformed state and $$\hat{\mathbf{{n}}}$$ is a unit vector normal to the tube wall. The expression (B.1) for the dimensional rate of working in the tube wall $${\mathrm{d}} E_{s}^{*}/ {\mathrm{d}} t^{*}$$ comes from integrating the product of the force from the transmural pressure and the normal component of the velocity of the tube wall, over the mid-plane of the tube wall. Inserting the appropriate limits for the integration within (B.1) and noting that $${p_{\mathrm{tm}}}^{*}$$ is independent of the azimuthal coordinate $$\tau$$ at leading order, it is found that   \begin{equation} \frac{{\mathrm{d}} E^{*}_{s}}{{\mathrm{d}} t^{*}} = \int_{0}^{L} {p_{\mathrm{tm}}}^{*} \int_{0}^{2\pi} \frac{\partial \mathbf{r}}{\partial t^{*}} \cdot \hat{\mathbf{{n}}} \, a h(\tau) \, {\mathrm{d}} \tau \, {\mathrm{d}} z^{*}. \label{eq:dimrow2} \end{equation} (B.2) Now, to leading order, the area perturbation can be written as   \begin{equation} A^{*} - A_{0}^{*} = \int_{0}^{2\pi} (\mathbf{r} - \mathbf{r}_{0}) \cdot \hat{\mathbf{{n}}} \, a h(\tau) \, {\mathrm{d}}\tau, \label{eq:areadefrel} \end{equation} (B.3) where $$\mathbf{r}_{0}$$ is the undeformed position of the tube wall mid-plane. By differentiating this relation with respect to $$t^{*}$$ and then substituting the result into (B.2), we obtain   \begin{equation} \frac{{\mathrm{d}} E^{*}_{s}}{{\mathrm{d}} t^{*}} = \int_{0}^{L} {p_{\mathrm{tm}}}^{*} \frac{\partial A^{*}}{\partial t^{*}} \, {\mathrm{d}} z^{*}. \label{eq:dimrow3} \end{equation} (B.4) The expression (B.4) is now non-dimensionalised using the scalings (2.10) for the axial length and cross-sectional area, (2.13) for the transmural pressure, (2.16) for the energy and $$t^{*} = Tt$$ for time. We find   \begin{equation} \frac{{\mathrm{d}} E_{s}}{{\mathrm{d}} t} = {\mathit{St}}^{2} \ell^{2} \int_{z_{1}}^{z_{2}} {p_{\mathrm{tm}}} \frac{\partial A}{\partial t} \, {\mathrm{d}} z, \label{eq:nondimrow1} \end{equation} (B.5) where $$E_{s}$$ is the dimensionless energy in the tube wall. We note that $$\partial A/\partial t = 0$$ in the rigid regions of the tube and thus we only need to take the integral within (B.5) over the region $$z_{1}