Global stability analysis of flexible channel flow with a hyperelastic wall

J. Fluid Mech. (2022), vol.934, A28, doi:10.1017/jfm.2021.1131 Global stability analysis of flexible channel flow withahyperelasticwall M.A. Herrada1, S. Blanco-Trejo1,J.Eggers 2and P.S. Stewart3,† 1E.S.I., Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Spain 2School of Mathematics, University of Bristol, Fry Building, Bristol BS8 1UG, UK 3School of Mathematics and Statistics, University of Glasgow, Mathematics and Statistics Building, University Place, Glasgow G12 8QQ, UK (Received 29 January 2021; revised 20 October 2021; accepted 12 December 2021) We consider the stability of flux-driven flow through a long planar rigid channel, where a segment of one wall is replaced by a pre-tensioned hyperelastic (neo-Hookean) solid of finite thickness and subject to a uniform external pressure. We construct the steady configuration of the nonlinear system using Newton’s method with spectral collocation and high-order finite differences. In agreement with previous studies, which use an asymptotically thin wall, we show that the thick-walled system always has at least one stable steady configuration, while for large Reynolds numbers the system exhibits three co-existing steady states for a range of external pressures. Two of these steady configurations are stable to non-oscillatory perturbations, one where the flexible wall is inflated (the upper branch) and one where the flexible wall is collapsed (the lower branch), connected by an unstable intermediate branch. We test the stability of these steady configurations to oscillatory perturbations using both a global eigensolver (constructed based on an analytical domain mapping technique) and also fully nonlinear simulations. We find that both the lower and upper branches of steady solutions can become unstable to self-excited oscillations, where the oscillating wall profile has two extrema. In the absence of wall inertia, increasing wall thickness partially stabilises the onset of oscillations, but the effect remains weak until the wall thickness becomes comparable to the width of the undeformed channel. However, with finite wall inertia and a relatively thick wall, higher-frequency modes of oscillation dominate the primary global instability for large Reynolds numbers. Key words: flow-vessel interactions †Email address for correspondence: [email protected] © The Author(s), 2022. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited. 934 A28-1 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart 1. Introduction Human physiology includes a wide number of examples of fluid flow through flexible-walled conduits including blood flow through the circulation (from rapid flow in the heart and large arteries to slow viscous flows through the capillaries), air flow through the lungs and upper airways, urine flows in the excretory system and peristaltic flows through the colon. In some circumstances these flows can exhibit instability, where the flow can interact with the flexible wall in a non-trivial way. Of particular interest in this study is the onset of self-excited oscillations, where the flow and the wall can spontaneously transition to an oscillatory limit cycle; in some cases this oscillation can even become chaotic. These oscillations manifest in physiological problems such as blood pressure measurement in the form of audible Korotkoff noises (Bertram, Raymond & Butcher 1989), and wheezing in the lung airways (Gavriely et al. 1989). Self-excited oscillations in flexible-walled vessels can be studied experimentally using a Starling resistor, a deceptively simple device featuring liquid flow driven through a section of externally pressurised flexible tubing mounted between two rigid pipes. Originally used as a flow resistor in cardiac experiments (Knowlton & Starling 1912), it has since become a canonical experiment for investigating fluid–structure interaction in its own right. In these experiments flow is driven using either a prescribed pressure or a prescribed flow rate, and the choice of set-up heavily influences the structure of the resulting oscillations. Results from the experiments are well summarised elsewhere (e.g. Bertram 2003; Grotberg & Jensen 2004; Heil & Hazel 2011), but we note that these self-excited oscillations occur in distinct frequency bands (Bertram, Raymond & Pedley 1990), and exhibit complicated nonlinear limit cycles which can be characterised using the methods of dynamical systems (Bertram, Raymond & Pedley 1991). Note that these experiments are typically conducted with relatively thick-walled tubes. For example, Bertram et al. (1990,1991) used tubes of wall thickness to baseline radius ratio of 0.3, while Bertram & Castles (1999) used tubes with a thickness to radius ratio of 0.37. There have been a number of theoretical studies of the Starling resistor set-up in an attempt to explain the underlying mechanisms leading to these different families of oscillation. Formulation of the full three-dimensional fluid structure interaction problem in a collapsible tube involves coupling unsteady Newtonian flow to a fully deformable elastic tube. While most theoretical models treat the tube wall as a thin shell, slightly reducing the complexity of the system, these models still require vast computational resources to resolve the unsteady oscillatory flow (Heil & Boyle 2010). Some analytical progress can be made in the limit of large membrane tension (where oscillations are high frequency, Whittaker et al. 2010), but this formulation is restricted to a state where the tube wall is almost uniform that has not yet been realised experimentally. The flexible tubing used in Starling resistor experiments is typically much thicker than is appropriate to model using thin shell theory. To date, the only theoretical studies which incorporate a thick-walled tube have been restricted to steady flow configurations (Marzo, Luo & Bertram 2005; Zhang, Luo & Cai 2018). In this paper we seek to address the stability of flow in a Starling resistor analogue with a thick hyperelastic wall, and investigate the role of wall thickness in promoting or inhibiting instability. Given the computational difficulty and expense of full three-dimensional unsteady models, theoretical study has often focused on empirical lumped parameter or cross-sectionally averaged models for flow in collapsible tubes (e.g. Shapiro 1977; Bertram&Pedley1982; Jensen 1990; Armitstead, Bertram & Jensen 1996), which have replicated many of the features noted in Starling resistor experiments, such as non-uniform steady profiles and spontaneous transition to self-excited oscillations in distinct 934 A28-2 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall oscillation frequencies. However, the flow field in these models is still approximate and misses many of the subtleties of flow separation and energy dissipation. To make progress in understanding the mechanisms of instability driving self-excited oscillations, a compromise system is needed which is less complicated than fully three-dimensional flow, but reduces the number of empirical assumptions needed for the lumped models. Pedley (1992) proposed a two-dimensional analogue of the Starling resistor, consisting of a planar rigid channel where a section of one wall has been replaced by a flexible sheet. This set-up has since become the subject of a wide variety of computational (e.g. Luo & Pedley 1995,1996,1998,2000;Heil2004) and theoretical studies (e.g. Jensen & Heil 2003; Guneratne & Pedley 2006;Stewartet al. 2010; Pihler-Puzovi´ c & Pedley 2013). Despite reduced computational cost compared with the three-dimensional tube system, a full exploration of the parameter space for this collapsible channel analogue has not yet been attempted, although progress toward quantifying the mechanisms of instability has been made in various regions of the parameter space. For example, in the case of prescribed upstream flux (the subject of this study), Xu, Billingham & Jensen (2014) quantified the mechanism driving ‘sawtooth’ oscillations in the asymptotic limit of a long downstream rigid section, where the nonlinear oscillation is driven by the resonance of two distinct modes of perturbation (mode-1 and mode-2) of similar frequency and the same wavelength, coupled by sloshing flow in the downstream rigid section. Furthermore, Huang (2001) simplified the flux-driven collapsible channel system by imposing an external pressure gradient on the flexible wall, which facilitated decomposition of the oscillatory flow into a sum of sinusoidal modes. This analysis reveals an alternative mechanism of oscillatory instability, driven by an imbalance between (unstable) downstream propagating waves (which transfer energy from the flow to the wall) and (stable) upstream propagating waves (which transfer energy back from the wall to the fluid). Further insights into the mechanisms of instability in these collapsible channel flows have been obtained using approximate one-dimensional models of the asymmetric channel system (derived using a flow-profile assumption, Stewart, Waters & Jensen 2009;Stewart et al. 2010; Xu, Billingham & Jensen 2013;Xuet al. 2014; Xu & Jensen 2015; Stewart 2017). In particular, a detailed exploration of the parameter space for flux-driven oscillations with constant external pressure was presented by Stewart (2017), where he found that when the fluid is inviscid, steady states only exist above a critical value of the membrane tension (for all other parameters held fixed), with a stable branch and an unstable branch (where the unstable branch is more collapsed than the stable branch). This critical point appears to be an organising centre of the dynamical system, in that many of the unsteady features of the system originate close to this point (such as the neutral curves for the two different families of self-excited oscillations). The importance of the critical point for inviscid steady states has previously been elucidated by Xu et al. (2013), who used an external pressure gradient. Stewart (2017) also described another branch of steady solutions maintained by viscous effects, which becomes increasingly collapsed as the wall tension is reduced. As the Reynolds number increases this viscous branch of steady solutions merges with one of the (essentially) inviscid branches. When the viscous branch merges smoothly with the stable inviscid branch then the stable steady state is unique. However, the other possibility is that the viscous branch merges with the unstable inviscid branch in a limit point bifurcation, where the system then exhibits three co-existing steady states across a narrow region of the parameter space: the stable inviscid solutions become the upper branch, the unstable inviscid solutions become the intermediate branch and the stable viscous solutions become the lower branch. 934 A28-3 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart Stewart (2017) also showed that the lower branch of steady solutions can become unstable to two distinct families of self-excited oscillation, with high and low frequency, respectively. However, in addition to the flow-profile assumption, this study considered the flexible wall to be a thin (massless) pre-stressed membrane with no bending rigidity. To overcome these simplifications, this study revisits the predictions of Stewart (2017)by modelling the flexible wall as a pre-tensioned hyperelastic solid, using the finite element method to compute the fully two-dimensional steady wall and flow profiles, and test their stability to time-dependent perturbations using a fully two-dimensional eigensolver. Our new model includes the wall thickness and wall mass as explicit parameters, and we investigate their influence on the predictions below. Another approach for theoretical modelling of this collapsible channel system has very recently been presented by Wang, Luo & Stewart (2021a,b), who treat the flexible wall as an asymptotically thin beam with resistance to both bending and stretching but with no pre-tension (based on an earlier model by Cai & Luo 2003;Luoet al. 2008). Using fully nonlinear simulations of this model, they identified a similar three-branch steady system for some parameters, showing that both the upper and lower branches of oscillation could (independently) become unstable to self-excited oscillations (Wang et al. 2021a) and these families of oscillations could merge together for low external pressures (Wang et al. 2021b). In this case the upper branch instability is restricted to a region in the near neighbourhood of that which exhibits multiple steady states (Wang et al. 2021b). In this study we also isolate a family of upper branch instabilities, but show that these are not limited to the region with multiple steady states but are instead unstable well away from the region of parameter space which exhibits instabilities of the lower steady branch (see §3.4 below). The role of wall mass in the onset of self-excited oscillations in flexible-walled vessels has already been considered for the flexible wall modelled as a thin membrane. For example, in the asymmetric channel system, Luo & Pedley (1998) coupled the heavy membrane to fully two-dimensional (unsteady) flow, showing that increasing the wall mass expands the region of parameter space where the system exhibits the primary global instability, and also results in an additional high-frequency oscillatory mode (superimposed on the fundamental mode) which eventually grows to dominate the lower-frequency mode. Also, Pihler-Puzovi´ c & Pedley (2014) investigated this channel system using interactive boundary layer theory, showing that wall mass drives an oscillatory instability which is always unstable in the presence of a cross-stream pressure gradient across the core flow (the system is always neutrally stable with no cross-stream gradient). Finally, Walters, Heil & Whittaker (2018) considered the role of wall mass in a thin shell model of flow in a collapsible tube in the limit of large pre-stress (where the tube is almost uniform), finding that wall inertia destabilises the primary mode of instability of the system while also lowering the corresponding oscillation frequency. In this paper we consider the planar channel analogue of the Starling resistor introduced by Pedley (1992), and propose a new numerical method to solve the combined fluid and solid problem based on that developed by Snoeijer et al. (2020) (which already has application to viscoelastic fluids, Eggers, Herrada & Snoeijer 2020). The model formulation is described in § 2, highlighting the novel features of the numerical method. In particular, we treat the elastic solid as a pre-tensioned hyperelastic material of uniform initial thickness with non-negligible density and subject to a uniform external pressure. We validate this numerical method against the steady predictions of Heil (2004), who considered an identical set-up with a thin shell model for the wall (§ 3.1), use unsteady simulations to examine the transition between the upper and lower branches of steady solutions (§ 3.2), examine the onset of self-excited oscillations from these steady solutions 934 A28-4 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall External gas Elastic solid Rigid solid h1 (x, t) L1L2 L p = 0 p = Pext h h2 (x, t) η1 ρ1 x y q e Fluid Ω1 (t) Ω2 (t) μ2 ρ2 Figure 1. Sketch of the flow geometry considered in this study. (§ 3.3), before using our new model to examine the role of membrane pre-tension (§ 3.4), the dynamics of oscillations growing from the upper branch of steady solutions (§3.5)as well as the role of wall thickness (§3.6) and wall inertia (§ 3.7) on the nonlinear steady solutions and the accompanying onset of oscillation. 2. Model formulation We consider the configuration sketched in figure 1, where an incompressible Newtonian fluid is flowing through a planar rigid (two-dimensional) channel of uniform internal width h. An interior section of length Lis removed from the upper wall of the channel and replaced by a pre-tensioned elastic solid of (initially) uniform thickness e, subject to a passive external gas at uniform pressure, Pext. This elastic wall can be deformed by the load of the external gas and by the fluid traction. The rigid sections upstream and downstream of the compliant segment are of length L1and L2, respectively. In this case the flow is driven by a prescribed upstream flux q, while the fluid pressure at the downstream end of the channel can be set to zero without loss of generality. The stability of this fluid–structure interaction problem has already been studied extensively using reduced models for the elastic wall (e.g. Luo & Pedley 1996; Jensen & Heil 2003;Luoet al. 2008;Stewart2017). In this work, we model the wall as a continuum hyperelastic solid of finite thickness, with no simplifications or reductions. Our formulation is based on first-order elasticity (elastic strain energy function dependent on the strain tensor), which places some restrictions on the boundary conditions that can be imposed. 2.1. Equations of motion The fluid domain Ω1is described by the planar coordinates x=xex+yey,wherex parametrises the lower wall of the channel, with x=0 at the intersection between the upstream rigid segment and the compliant segment, while yparametrises the direction normal to the entirely rigid wall pointing into the fluid (in the plane of the channel). The solid domain Ω2is measured relative to a reference configuration parametrised by the 934 A28-5 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart coordinates X=Xex+Yey,whereXparametrises the lower surface of the flat wall and Y parametrises the direction pointing into the wall (in the plane of the channel). The conservation of mass and momentum equations in the fluid (i=1) and solid (i=2) subdomains are given by ∇·vi=0,(i=1,2), (2.1a) ρi∂vi ∂t+(vi·∇ )vi=∇·σi,(i=1,2), (2.1b) where ρiis the density, Iis the identity tensor, vithe velocity field and σiis the stress tensor of material i(i=1,2). Each stress tensor depends on the characteristics of the material through a constitutive model. In region 1 we consider an incompressible Newtonian fluid, where this stress tensor takes the form σ1=−p1I+η1∇v1+∇vT 1,(2.1c) where p1is the fluid pressure and η1is the fluid viscosity. In region 2 we consider a neo-Hookean (hyperelastic) solid which has a pre-stress, σ(0) 2p, in the initial undeformed state, where the stress tensor is given by (Snoeijer et al. 2020) σ2=−p2I+μ2F·FT−I+F·σ(0) 2p·FT,(2.1d) where p2is the solid pressure, μ2is the elastic shear modulus, x(X,t)is the position of a material point after deformation of the solid and F=∂x/∂Xis the deformation gradient tensor. In the initial state, x=Xand F·FT=I. To make a connection between the Eulerian formulation for the conservation of mass and momentum equations for the solid ((2.1)withi=2) and the Lagrangian formulation for the elastic stress, we need to determine the deformation generated by transport by the solid velocity v2. This is achieved using the inverse Lagrangian map X(x,t)(Kamrin, Rycroft & Nave 2012), which satisfies ∂X ∂t+v2·∇X=0,(2.1e) because the reference coordinates are invariant under the flow. Given the bi-dimensionality of the problem, the material points can be expressed in Cartesian coordinates and so the velocity vectors can be written as vi=vyiey+vxiex,(i=1,2), (2.1f) while the stress tensors can be written as σ=σyyey⊗ey+σyxey⊗ex+σxyex⊗ey+σxxex⊗ex,(2.1g) and finally the deformation tensor in the solid can be written as F=∂y ∂Y ey⊗ey+∂y ∂X ey⊗ex+∂x ∂Y ex⊗ey+∂x ∂X ex⊗ex.(2.1h) In the undeformed position the elastic solid is subject to an initial longitudinal tension, To, and therefore the initial stress is σ(0) 2p=(T0/e)ex⊗ex. 934 A28-6 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall For the elastic domain, it is convenient to replace the incompressibility equation based on the velocity field ((2.1a)withi=2) by a constraint involving the deformation tensor F (Snoeijer et al. 2020)intheform det(F)=∂y ∂Y ∂x ∂X−∂y ∂X ∂x ∂Y=1.(2.1i) To impose the upstream flux boundary condition for the liquid, we impose a Poiseuille profile at the channel entrance, x=−L1,intheform v1x=6q h3y(h−y), v1y=0,(x=−L1,0⩽y⩽h). (2.1j) At the channel exit, x=L+L2, we impose zero fluid pressure, p1=0. Along the entirely rigid wall we apply no-slip conditions in the form vx1=vy1=0,(y=0,−L1⩽x⩽L+L2). (2.1k) Similarly, along the rigid parts of the upper wall we apply no-slip boundary conditions in the form vx1=vy1=0,(y=h,−L1⩽x⩽0,x⩾L). (2.1l) We assume that the flexible surface (where the elastic solid and the fluid interact) can be written as a function of x(i.e. the surface does not overturn or expand beyond the range 0⩽x⩽L), so that y=h1(x,t). Across this interface we impose that the velocity field must be continuous, in the form vx1=vx2,v y1=vy2,(y=h1,0⩽x⩽L), (2.1m) and impose a balance of normal and tangential stresses between the solid and the fluid, in the form n1·(σ1−σ2)·n1=0,t1·(σ1−σ2)·n1=0,(2.1n) where n1=ey−exh1,x (1+h2 1,x)1/2,t1=ex+eyh1,x (1+h2 1,x)1/2,(2.1o) are normal and tangential vectors to the surface y=h1(x,t), respectively, and the subscript xrepresents a derivative with respect to x. In this first-order elasticity approach we enforce no deformation along the surfaces where the elastic material is adhered to the rigid walls (i.e. the displacement of the solid is clamped along two edges of the rectangle in contact with the rigid walls), in the form v2x=v2y=0,Y=y,X=x,(x=0,x=Lwith h⩽y⩽h+e). (2.1p) However, our approach does not replicate the resistance to bending of a classical Euler–Bernoulli beam. This would require second-order (or strain gradient) elasticity, where the elastic strain energy function is assumed to depend on both the strain tensor and the strain gradient tensor (Bertram & Forest 2020). In that case one must impose additional constraints on the contact between the beam and the rigid wall e.g. conditions on the derivatives of displacement, such as prescribed slope or torque. Finally, we denote 934 A28-7 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart the external surface of the flexible wall as y=h2(x,t),(0⩽x⩽L) and impose that the normal and tangential elastic stresses are balanced with the external pressure, in the form n2·(σ2−PextI)·n2=0,t2·(σ2)·n2=0,(2.1q) where n2=ey−exh2,x (1+h2 2,x)1/2,t2=ex+eyh2,x (1+h2 2,x)1/2,(2.1r) are normal and tangential vectors to the surface y=h2(x,t). 2.2. Mapping technique The numerical technique used in this study is a variation of that developed by Herrada & Montanero (2016) for interfacial flows and extended by Snoeijer et al. (2020)to apply to hyperelastic solids. The spatial domain occupied by the fluid, Ω1(t), is mapped onto a rectangular domain (parametrised by Cartesian coordinates ξ1and χ1,whereξ1 parametrises the lower rigid wall and χ1parametrises the channel inlet) by means of a non-singular mapping y=f1(ξ1,χ 1,t), x=g1(ξ1,χ 1,t), [−L1⩽ξ1⩽L+L2]×[0 ⩽χ1⩽1],(2.2) where the shape functions f1and g1are obtained as part of the solution. In order to capture large anisotropic deformations, the following quasi-elliptic transformation (Dimakopoulos & Tsamopoulos 2003) was applied: g22 ∂2f1 ∂ξ2 1 +g11 ∂2f1 ∂χ2 1 −2g12 ∂2f1 ∂ξ1∂χ1 =Q,(2.3a) g22 ∂2g1 ∂ξ2 1 +g11 ∂2g1 ∂χ2 1 −2g12 ∂2g1 ∂ξ1∂χ1 =0,(2.3b) where the coefficients take the form g11 =∂g1 ∂ξ12 +∂f1 ∂ξ12 ,g22 =∂g1 ∂χ12 +∂f1 ∂χ12 ,g12 =∂g1 ∂χ1 ∂g1 ∂ξ1 +∂f1 ∂χ1 ∂f1 ∂ξ1, (2.4a–c) with Q=−∂D1 ∂χ1 ∂f1 ∂ξ1 −∂D1 ∂ξ1 ∂f1 ∂χ1J D1,J=∂g1 ∂χ1 ∂f1 ∂ξ1 −∂g1 ∂ξ1 ∂f1 ∂χ1,(2.5a,b) and D1=p   ∂f1 ∂ξ12 +∂g1 ∂ξ12∂f1 ∂χ12 +∂g1 ∂χ12+(1−p). (2.6) In the above expressions, pis a free parameter between 0 and 1 where the case p=0 corresponds to the classical elliptical transformation. All the simulations in this work were conducted using p=0.2. Although there is no overturning in the wall profiles for the cases analysed in this work, this transformation of the liquid domain facilitates the analysis of more complicated geometries. For example, it has been successfully used to describe pinch-off in pendant drops (Ponce-Torres et al. 2020). 934 A28-8 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall The spatial domain occupied by the elastic solid in the current stage, Ω2(t), and in the initial stage, Ω2o, are also mapped onto rectangular domains (parametrised by Cartesian coordinates ξ2and χ2,whereξ2parametrises the lower surface of the flexible wall and χ2parametrises the edges in contact with the rigid segments of the channel) by means of non-singular mappings in the form y=f2(ξ2,χ 2,t), x=g2(ξ2,χ 2,t), Y=F2(ξ2,χ 2,t), X=G2(ξ2,χ 2,t), [0 ⩽ξ2⩽L]×[0 ⩽χ2⩽1],(2.7) where again the functions f2,g2,F2and G2should be obtained as a part of the solution. To determine these functions, the following equations have been used: g2=ξ2,(2.8a) F2=h+eχ2.(2.8b) Note that (2.8a) guarantees that the discretisation used for the variable ξ2is automatically applied to variable x. Finally, (2.8b) indicates that at the initial stage the elastic part of the upper channel wall is a perfect rectangle of uniform width e. Some additional boundary conditions for the shape functions are needed to close the problem. At the channel entrance, we impose g1=−L1,f1=hχ1,(x=ξ1=−L1), (2.9a) while at the channel exit, we use g1=L+L2,f1=hχ1,(x=ξ1=L+L2). (2.9b) On the lower wall, we impose g1=ξ1,f1=0,(y=χ1=0), (2.9c) while on the rigid parts of the upper channel wall, we use g1=ξ1,f1=h,(−L1⩽x=ξ1⩽0,x=ξ1⩾L,y=h). (2.9d) At the flexible surface, we also impose f1=f2,g1=g2,(0⩽x=ξ1=ξ2⩽L,y=h1(x,t), χ1=1,χ 2=0). (2.9e) Finally, we enforce no displacement of the elastic solid along the two edges of the rectangle in contact with the rigid walls, in the form g2=G2=ξ2,f2=F2=h+eχ2, (x=ξ2=0,x=ξ2=L,h⩽y⩽(h+e), 0⩽χ2⩽1).(2.9f) Figure 2 shows an example of the mappings used in this work. The green (magenta) lines represent the liquid (solid) mesh in the real space (bottom panel) and in the computational domain (top panel). The unknown variables in the liquid domain are f1,g1,p1,v1xand v1y while the unknown variables in the solid domain are f2,g2,p2,v2x,v2y,F2and G2.All the derivatives appearing in the governing equations are expressed in terms of χ,ξand t. These mappings are applied to the governing equations (2.1) and the resulting equations are discretised in the χ-direction with nχ1and nχ2Chebyshev spectral collocation points in the liquid and solid domains, respectively. Conversely, in the ξ-direction we use fourth-order finite differences with nξ1and nξ2equally spaced points in the liquid and solid domains, respectively. The results presented in this work were carried out using nξ1=641, nξ2=201, nχ1=19 and nχ2=14. In the Appendix we demonstrate that the eigenvalues characterising the linear modes do not change significantly when the number of grid points is increased. 934 A28-9 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart but as the channel becomes increasingly constricted the rate of collapse increases and boundary layer separation takes place (figure 6c), where a re-circulation region becomes evident close to the downstream outlet of the compliant segment of the channel (figure 6d), creating a region of much lower pressure (figure 6e). A movie showing the entire transition is provided in the online supplementary material available at https://doi.org/10.1017/jfm. 2021.1131. There is an interesting analogy between these observations and those reported for swirling flows in pipes (see for e.g. Lopez 1994; Herrada, Pérez-Saborid & Barrero 2003), where fluid flows with a significant azimuthal velocity component through a rigid circular tube with an axisymmetric (fixed) sinusoidal indentation over a finite length. In this analogy the indentation of the pipe mirrors the collapse of the compliant segment of the channel, while the azimuthal fluid velocity component (and to some extent the compressibility of the fluid) extracts energy from the mean flow in a similar way to the compliance of the elastic wall. These swirling flows exhibit multiple (stable) steady solutions for a given set of parameters (when the Reynolds number is larger than a critical one) and the steady solutions can be described using bifurcation diagrams with three branches of steady solutions and two limit points, analogous to those presented in figure 4; this behaviour was recently termed ‘double hysteresis’ (Shtern 2018). These swirling flows also exhibit an unsteady transition from a nearly columnar flow to a recirculating flow when the swirling parameter is larger than a critical value (vortex breakdown), analogous to the spontaneous collapse of the channel we observe as the external pressure increases above the critical value (ˆpext1). In the former case, centrifugal forces generate an adverse axial pressure gradient that induces a recirculating flow, whereas the channel collapse generates an adverse pressure gradient that drives detachment of the boundary layer adjacent to the flexible wall. The flow structures in figures 6–9 of Herrada et al. (2003) are reminiscent of the transition observed in figure 6, where in both cases the vortex breakdown occurs just downstream of the point of greatest indentation. The only significant difference comes in the cross-stream location of vortex shedding: the symmetry of the cylindrical geometry in the swirling flows results in vortex shedding near the axis of the tube, while in the collapsible channel the vortex shedding occurs near the flexible wall. 3.3. Linear stability results Having computed the steady configurations of the system, we now analyse the temporal linear stability of the three different steady solution branches depicted in figure 4. For this large value of pre-tension (ˆ T0=10) we find that the steady solutions along the section of the upper branch tested are globally stable to time-dependent perturbations (all the eigenvalues have ωi<0) for all external pressures greater than the outlet pressure (i.e. ˆpext ⩾0), while the solutions along the intermediate branch are always unstable (at least one eigenvalue has ωi>0withωr=0). Figure 7 illustrates the stability of the lower steady branch, showing the eigenvalue spectrum of the frequency ωfor several values of the external pressure, ˆpext. In this case (and in figure 11 below) we focus only on the most unstable eigenvalues, illustrating those with ωi>−0.5.We find that the lower branch is stable for sufficiently small external pressure, becoming globally unstable via a Hopf bifurcation when the external pressure exceeds a critical value, ˆp∗ ext ≈1.752 (i.e. a pair of complex conjugate eigenvalues cross the real axis with non-zero ωr). At this critical point, the corresponding steady state is shown in figure 7(b), where it is inflated at the upstream end and collapsed at the downstream end (termed mode-2). The corresponding eigenfunction of the wall profile for the neutrally stable mode is shown in figure 7(c), which has two extrema (mode-2). We label the oscillatory modes associated with the 934 A28-16 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall (a)(b) (c) h ˆ 1b ωˆi  δh1 00 0.5 1.0 501015 –1.0 51015 0 1.0 0.5 –0.5 –2 –1 0 1 2 –0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 mode-(i), Re(δh1) mode-(i), Im(δh1) pˆext = 1.6 pˆext = 1.752 pˆext = 2 pˆext = 2.2 pˆext = 2.5 ωˆr xˆ   Figure 7. Stability of the lower steady branch to time-dependent perturbations for fixed Reynolds number (Re =500) and fixed pre-tension (ˆ T0=10): (a) five eigenvalue spectra in the ω-plane for increasing values of ˆpext;(b) profile of the lower surface of the steady wall at neutral stability (ˆpext ≈1.752); (c) real and imaginary parts of the wall profile eigenfunction at neutral stability (ˆpext ≈1.752). Here,ˆe=0.01 and ˆρ=0. lower branch with lower case Roman numerals (i), (ii), (iii) ... in the order of increasing frequency, which is generally the order they become unstable as the external pressure increases, and so this primary instability is denoted mode-(i). These stability predictions agree well with the results presented by Heil (2004), where his figure 5 shows that the flow becomes unsteady and exhibits self-excited oscillations for ˆpext =2.5, well inside our unstable regime. These results are also qualitatively similar to the predictions of the one-dimensional model of Stewart (2017), who showed that his lower branch of steady solutions becomes unstable to a mode-2 oscillation as the primary global instability of the system as the external pressure increases. We overview the parameter space in figure 8 to summarise the regions of interest. We illustrate the region with multiple steady solutions by tracing the value of the external pressure at the limit points of the upper and lower steady branches (ˆpext1and ˆpext2, analogous to those found in figure 4) as a function of the Reynolds number; similar to Stewart (2017), we find that this region with multiple steady states exists for Reynolds numbers greater than a threshold (Re >Recusp ≈330). We further plot the critical external pressure for the onset of oscillatory instability, ˆp∗ ext, as a function of the Reynolds number, finding that for the range of Reynolds numbers explored here the neutral stability curve lies entirely within the range where there is a unique steady solution along the lower steady branch, so ˆp∗ ext >ˆpext1. Note that we observe no instability of the upper steady branch for this choice of the wall pre-tension (ˆ T0=10) across the range 0 ⩽ˆpext ⩽ˆpext1.Itemerges below that this branch only becomes unstable for ˆpext <0 for this value of ˆ T0, which is not considered here. For large Reynolds number we might expect the neutral stability curve to enter the region of parameter space with multiple steady states (in a similar manner to Stewart 2017), but this possibility is discussed in more detail below. 3.4. The influence of the pre-tension in the solid When the pre-tension of the elastic wall is reduced, we observe a decrease in the critical Reynolds number beyond which multiple steady flows exist, and the steady state 934 A28-17 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart p ˆext Re Recusp ≈ 330 50 100 150 200 250 300 350 400 Multiple steady flows 500450 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 Self-excited oscillations Stable pˆext2 pˆext1 mode-(i), pˆ∗ ext Figure 8. Overview of the critical conditions for self-excited oscillations for pre-tension ˆ T0=10, plotting the critical external pressure for instability as a function of the Reynolds number. The cross symbol indicates the point in parameter space which corresponds to the unsteady simulation shown in figure 6.Here ,ˆe=0.01, ˆρ=0. bifurcation diagram and neutral stability curves become more complicated. To illustrate this complexity, in figure 9 we characterise the multiplicity of steady solutions that exist for a lower value of the pre-tension (ˆ T0=5) while holding the Reynolds number fixed (Re =500), plotting the minimal (ˆ hmin) and maximal (ˆ hmax) widths of the steady channel as a function of the external pressure, for the upper and lower branches of steady solutions, obtained following the same procedure as §3.1. Similar to the case we considered in figure 4 (ˆ T0=10), when the external pressure increases beyond a certain value, ˆpext =ˆpext1, there is a jump from a solution on the upper branch to a solution on the lower branch (where the channel becomes much more collapsed). In the same way, as we decrease the external pressure along the lower branch below a certain value, ˆpext =ˆpext2, there is a jump back to the upper branch. To overview these steady solutions across the parameter space, in figure 10 we plot the external pressure at the limit points of the steady solutions (ˆpext1and ˆpext2) as a function of the Reynolds number for a lower value of the pre-tension (ˆ T0=5), where we find that the critical Reynolds number for multi-valued solutions has reduced (Recusp ≈275inthis case). To further illustrate the stability of these steady solutions, in figure 10 we also trace the critical external pressure for the onset of instability as a function of the Reynolds number, finding again that the lower branch of steady solutions (branch III) becomes unstable for external pressures greater than a critical value, ˆp∗ ext, and is stable otherwise (figure 10). This observation is similar to our observation for large pre-tension (ˆ T0=10), with the only difference that now the loss of stability is closer to the region of multiplicity of steady solutions, with the two bounding curves almost overlapping for the largest Reynolds numbers considered. Tracing these curves to larger Reynolds numbers is an interesting direction of future work, where we might expect the neutral stability curve and the trace of the lower branch limit point to eventually intersect. Such an intersection was previously observed by Stewart (2017), where the Hopf bifurcation (associated with the 934 A28-18 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall pˆ∗ extI 0.5 1.0 1.5 2.0 h ˆ max h ˆ min pˆext2pˆext1pˆ∗ ext 1.2 1.3 1.4 1.5 1.6 Stable Self-excited oscillations Self-excited oscillations 1.7 1.8 h ˆ max h ˆ min (I) (I) h ˆ max h ˆ min (II ) (II ) h ˆ max h ˆ min (III ) (III ) pˆext Figure 9. Nonlinear steady solutions of the model for fixed Reynolds number (Re =500) and pre-tension (ˆ T0=5), showing the maximal and minimal channel widths as a function of the external pressure ˆpe.Here, ˆe=0.01 and ˆρ=0. Self-excited oscillations Self-excited oscillations + Multiple steady flows p ˆext Re Recusp ≈ 275 100 200 300 400 500 600 700 800 1.0 0.5 0 1.5 2.0 2.5 3.0 Stable pˆext2 pˆext1 mode-(i), pˆ∗ ext mode-(a), pˆ∗ extI Figure 10. Overview of the critical conditions for self-excited oscillations for lower pre-tension ˆ T0=5, plotting the critical external pressure for instability as a function of the Reynolds number. The plus symbol indicates the point in parameter space which corresponds to the nonlinear portrait of the upper branch instability shown in figure 12.Here ,ˆe=0.01, ˆρ=0. oscillation) and the saddle node bifurcation (associated with the steady solutions) interact in a co-dimension 2 bifurcation, suggesting a nearby homoclinic orbit (Glendinning 1994). However, for this lower value of the pre-tension we also observe that steady solutions along the upper branch (branch I in figure 9) also become temporally unstable for external pressures below a critical value, denoted ˆp∗ extI, and are stable otherwise (see figure 10). This means that for Re >Recusp there is only a narrow interval of external pressures compatible with a steady stable flow, focused around the region with multiple steady solutions. Instability of the upper branch of steady solutions has recently been noted by 934 A28-19 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart ωˆr pˆext = 0.7 pˆext = 1.12 pˆext = 1.3 pˆext = 1.5 (a)(b) (c) h ˆ 1b ωˆi  δh1 00 0.5 1.5 1.0 50 10 15 –1.0 51015 0 1.0 0.5 –0.5 –1.0 –0.5 0 0.5 1.0 –0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 mode-(a), Re(  δh1) mode-(a), Im(  δh1) xˆ Figure 11. Exploration of the lower branch instability for lower pre-tension ˆ T0=5: (a) five eigenvalue spectra in the ωplane for increasing external pressure; (b) steady wall profiles for the choice of external pressure where the system is neutrally stable; (c) real and imaginary parts of the corresponding eigenfunction of the wall profile at neutral stability. Here,ˆe=0.01 and ˆρ=0. Wang et al. (2021a) using a flexible wall modelled as a thin nonlinear beam, but in their case the region of instability is located within and directly adjacent to the region with multiple steady solutions (Wang et al. 2021b), in contrast to that noted here. The fully developed limit cycles also exhibit some significant differences (see § 3.5 below). To further explore this upper branch instability for lower pre-tension (ˆ T0=5) and fixed Reynolds number (Re =400), in figure 11(a) we plot the corresponding eigenvalue spectra for several values of the external pressure, where a complex conjugate pair of eigenvalues cross into the upper half-plane for ˆp<ˆp∗ extI ≈1.12 (Note that ˆp∗ ext ≈1.66), consistent with a Hopf bifurcation. At neutral stability the steady configuration of the flexible wall is entirely inflated with a single hump (termed mode-1, see figure 11b), while the neutrally stable eigenfunction of the oscillating wall profile is mode-2 (figure 11c), similar to the instability of the lower branch. Note that the frequency of oscillation along the upper branch is generally larger than the corresponding instability along the lower branch. Given that this oscillation also has a mode-2 structure of the wall shape eigenfunction, we label modes associated with the upper branch using Roman letters (a),(b),... in order of increasing frequency, which is generally the order they become unstable as the Reynolds number increases. The primary oscillatory mode associated with the upper branch is therefore labelled mode-(a). It is interesting to note that the instability of the mode-1 steady state exhibits a mode-2 eigenfunction profile, presumably because the prescribed upstream flux suppresses modes that induce large volume changes in the flexible segment of the channel (such as the mode-1 oscillations observed with prescribed upstream pressure e.g. Jensen & Heil 2003;Stewartet al. 2009,2010). The upper branch neutral stability point, ˆp∗ extI, can be traced (by numerical continuation) to larger values of the wall pre-tension; we find that the critical external pressure must become negative to induce instability for ˆ T0=10, explaining why it was not observed in figures 7 and 8, where we restrict attention to external pressures larger than the channel outlet pressure (ˆpext >0). 934 A28-20 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall We note that the neutral stability curves associated with both the upper and lower steady branches trace close to the region with multiple steady solutions as the Reynolds number increases, suggesting this region plays a key role in the structure of the dynamical system. Stewart (2017) showed that the limit point on the upper steady branch (traced by the blue curve in figures 8 and 10) asymptotes to the saddle node bifurcation point for steady solutions of the inviscid system as the Reynolds number increases. Indeed, both Xu et al. (2013)andStewart(2017) identified the threshold where inviscid steady states emerge as an organising centre of the dynamical system, consistent with our observation. Conversely, the lower branch of steady solutions is entirely maintained by the fluid viscosity (Stewart 2017), and is thus absent in the inviscid limit. 3.5. Limit cycles of upper branch instability Fully nonlinear simulations of self-excited oscillations growing from the lower branch of steady solutions have been widely reported elsewhere (e.g. Heil 2004;Luoet al. 2008). An instability of the upper branch of steady solutions was recently reported by Wang et al. (2021a), who considered flow through a similar two-dimensional collapsible channel system modelling the flexible wall as a thin (nonlinear) beam with resistance to both bending and stretching (with no pre-stress), and the nonlinear limit cycles were explored using fully nonlinear simulations. However, the upper branch oscillations evident from the present model exhibit a significant difference in structure: for the oscillations reported by Wang the unstable region restabilises as the upper branch limit point is reached (Wang et al. 2021a) and remains confined to the neighbourhood of the region with multiple steady states (Wang et al. 2021b), whereas for the present model the system is stable in the neighbourhood of the upper branch limit point and instead the unstable region extends over a wide range of external pressures away from the region with multiple steady states (figure 10). Given the difference in structure between our predictions and those of Wang et al. (2021a), in figure 12 we examine the underlying dynamics of our upper branch oscillations using fully nonlinear simulations of our model (method described in § 2.5) at a point in parameter space within the upper branch neutral stability curve. In this case we choose Re =500, ˆpext =1andˆ T0=5, marked with a plus inside the unstable region in figure 10. Initiating the simulation on the upper branch steady solution, numerical noise is enough to trigger an oscillatory instability evident in the time trace of the maximal channel width (see figure 12(a) with growth rate and frequency consistent with the global linear stability eigensolver), eventually saturating into a complicated nonlinear limit cycle (one period shown in figure 12b). A movie showing the flow field and vorticity over several periods of this limit cycle is provided in the online supplementary material. Over a period of this limit cycle the wall profile grows a single hump at the downstream end of the compliant segment (figure 12c); this hump propagates upstream reaching a global maximum (figure 12d) before being reflected back downstream again by the upstream rigid segment, where its amplitude subsequently decreases. As this hump propagates downstream a second hump appears at the downstream end of the compliant segment (figure 12e) which eventually dominates the first (figure 12f). However, these two humps do not coalesce but instead the x-location of the maximum wall deflection changes discontinuously at the global minimum of ˆ hmax (figure 12(f), explaining the cusp in figure 12(b)att≈1138.2,1145.5,1152.8). This second hump grows in amplitude, engulfing the remains of the first hump and shedding a low pressure vortex into the downstream rigid segment (figure 12g). This propagating vortex creates a so-called 934 A28-21 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart (c)(e) (g)(h) (d) (f) 1.6 1.8 2.0 2.2 (a) (b) h ˆ max h ˆ max 0 200 400 800600 Linear regime Fully developed nonlinear regime 1000 1138 1140 1142 1144 1146 1148 t 1150 1152 1154 1156 1158 1.8 1.9 2.0 2 1 0 (g) 500 –500 0 500 –500 0 500 –500 0 pˆbpˆb 0510 t = 1155 15 2 1 0 (h) 500 –500 0 0510 t = 1155.4 15 2 1 0 (e) 0510 t = 1152.8 15 2 1 0 (f) 500 –500 0 0510 t = 1154 15 2 1 0 (c) 0510 t = 1149.6 15 2 1 0 (d) 500 –500 0 0510 t = 1151.2 15 xˆxˆ yˆyˆ yˆyˆ yˆyˆ Figure 12. The mechanism of upper branch instability for a thin hyperelastic wall (ˆe=0.01) with no wall inertia ( ˆρ=0): (a) the maximal channel width ˆ hmax as a function of time; (b) zoom-in over panel (a) over one period of oscillation; streamlines and pressure colour map of the channel close to the outlet of the compliant segment at six selected times over a period of oscillation including: (c)t=1149.6; (d)t=1151.2; (e)t= 1152.8; ( f)t=1154; (g)t=1155; (h)t=1155.4. The fully developed limit cycle of interest is enclosed in the red box in (a). The times corresponding to the snapshots in panels (c–h) are labelled in (b). Here,Re =500, ˆpext =1andˆ T0=5. vorticity wave in the downstream rigid segment (particularly evident in figure 12c,g,h) while the large hump at the downstream end of the compliant segment drives a short region of channel collapse at the upstream end. As this vorticity wave propagates downstream, the single hump in the compliant segment propagates upstream, repeating the cycle. The nature of this oscillation exhibits many of the features of the nonlinear upper branch oscillations described by Wang et al. (2021a), including the development of an upstream propagating hump. However, for their upper branch oscillations this hump is annihilated by the upstream rigid segment (not reflected) and the flow remains entirely laminar throughout, with no evidence of low pressure vortex shedding. However, the present model is restricted by the assumption of first-order elasticity, meaning that we cannot apply as 934 A28-22 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall h ˆ min νˆmax 0 –0.10 0.90 0.95 1.00 –0.05 (c) 0 0.5 1.0 2.01.5 0 0.5 1.0 2.01.5 0 0.5 1.0 2.01.5 0 0.5 1.0 2.01.5 Self-excited oscillations Stable ωˆiωˆr 0.50 0.38 0.40 0.42 0.48 (a) (d) (b) eˆeˆ Figure 13. The influence of the wall thickness on the steady and oscillatory solutions in the absence of wall inertia ( ˆρ=0): (a) the minimal steady channel width ˆ hmin as a function of the wall thickness; (b) the maximal steady flow speed vmax as a function of the wall thickness; (c) the growth rate of the primary oscillatory mode as a function of the wall thickness (mode-(i)); (d) the frequency of the primary oscillatory mode (mode-(i)) as a function of the wall thickness. Here,ˆ T0=5, ˆpext =2.98 and Re =50. many boundary conditions at each end of the beam as Wang et al. (2021a) (who applied zero slope conditions at the end of the beam in addition to the clamped conditions). These vorticity waves have previously been observed in channel flows with self-excited oscillations from a collapsed (lower branch) steady state (Luo & Pedley 1996;Luoet al. 2008) or with prescribed (oscillatory) wall motion in one compartment (Stephanoff et al. 1983; Pedley & Stephanoff 1985). 3.6. The influence of the wall thickness In this subsection we analyse the influence of the dimensionless wall thickness, ˆe,onthe model predictions. We consider a particular case holding the pre-tension, external pressure and Reynolds number fixed (ˆ T0=5, ˆpext =2.98 and Re =50). For these parameters, with wall thickness ˆe=0.01, the system has a unique steady wall shape where the external pressure is sufficiently large to collapse the channel wall (ˆ hmin <1). These parameters are chosen so that the system is just inside the unstable regime for lower branch oscillations (Re =50 and ˆ T0=5 which has critical ˆp∗ ext ≈3.001). In figure 13 we characterise how an increase in the wall thickness influences the underlying steady flow (figure 13a,b)andthe critical conditions for the onset of lower branch oscillations (figure 13c,d). Considering the steady system first, figure 13(a) shows that the increase in wall thickness has little effect on the overall shape of the flexible wall for this value of Reynolds number; the channel becomes slightly less constricted as the wall thickness increases. Similarly, figure 13(b) shows that increasing wall thickness slightly reduces the maximal streamwise velocity through the constriction, ˆvmax =maxx,y(ˆv1xb)(as expected by conservation of mass). However, the wall thickness plays a more significant role in determining the stability of these steady solutions. The increase of the wall thickness results in the initially unstable solution (for ˆe=0.01) becoming stable for a critical value of the wall thickness ˆe0.08 (figure 13c), with a corresponding decrease in the frequency of oscillation (figure 13d). In order to quantify the effect of increasing the wall thickness on the stability of the system across the parameter space, in figure 14 we plot the critical external pressure 934 A28-23 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart ωˆr 50 100 150 200 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 mode-(i), eˆ = 2 mode-(i), eˆ = 1 mode-(i), eˆ = 0.4 mode-(i), eˆ = 0.2 mode-(i), eˆ = 0.01 mode-(i), eˆ = 2 mode-(i), eˆ = 1 mode-(i), eˆ = 0.4 mode-(i), eˆ = 0.2 mode-(i), eˆ = 0.01 p ˆ∗ ext Self-excited oscillations Stable Re 50 100 150 200 Re (a)(b) Figure 14. Overview of the critical conditions required for the onset of oscillations while changing the wall thickness, plotting the critical external pressure for the onset of oscillation ˆp∗ ext as a function of the Reynolds number for five different wall thicknesses (ˆe=0.01, ˆe=0.2, ˆe=0.4, ˆe=1andˆe=2). Here,ˆ T0=5. for the onset of the primary oscillatory instability of the lower branch (mode-(i)), denoted as ˆp∗ ext, as function of the Reynolds number for fixed pre-tension (ˆ T0=5) and five different wall thicknesses (ˆe=0.01,0.2,0.4,1,2). Note that we have limited our investigation to values of the Reynolds numbers smaller than the critical value required for multiple steady solutions (Recusp), so the steady profile is unique. In the absence of wall inertia (ˆρ=0)the effect of the wall thickness on the critical conditions for instability remains weak for relatively thin walls (ˆe=0.01,0.2,0.4): the steady flow remains almost unchanged (figure 13a,b) and there is only a mild stabilisation of the instability, characterised by an increase in the critical pressure needed to generate self-excited oscillations (figure 14). It emerges that the thickness of the wall must be of the order of the channel width (i.e. ˆe∼1) before there is any significant difference in the stability threshold. For example, for ˆe=1andˆe=2the critical pressure for the onset of instability is more appreciably increased compared with ˆe=0.01 (figure 14a), while the oscillation frequency is decreased (figure 14b). Furthermore, for ˆe=2 the critical external pressure and oscillation frequency both saturate as the Reynolds number becomes large (figure 14). We show in § 3.7 below that changes to the stability of the system are even more prominent when we include wall inertia. 3.7. The influence of wall inertia We now examine the influence of increasing wall inertia. It should be noted that the steady version of the full nonlinear equations (2.1) is independent of the wall inertia parameter ˆρ, and so all steady results are unchanged from those reported above. To study the additional influence of wall inertia on the onset of self-excited oscillations growing from the lower branch of static solutions, in figure 15 we trace the growth rate (figure 15a) and frequency (figure 15b) of the mode-(i) instability from figures 8 934 A28-24 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press Flexible channel flow with a hyperelastic wall ωˆi ωˆr (a) (b) (c) (d) (e) –1 501015 5 010 15 501015 –1 –1 0 1 0 1 0 1 2 3 4 0 1 0 20406080100 20 40 60 80 100 –0.10 –0.05 –0.15 0 0.15 0.10 0.05  δh1  δh1  δh1 mode-(i) mode-(ii) mode-(iii) mode-(iv) mode-(i) mode-(ii) mode-(iii) mode-(iv) mode-(i), Re(  δh1) mode-(i), Im(  δh1) mode-(ii), Re(  δh1) mode-(ii), Im(  δh1) mode-(iii), Re(  δh1) mode-(iii), Im(  δh1) ρˆxˆ Figure 15. The role of increasing wall inertia in the growth rate and frequency of self-excited oscillations for fixed wall thickness ˆe=0.2: (a) the growth rate of the first four oscillatory modes as a function of the wall inertia parameter; (b) the corresponding frequency of the first four oscillatory modes as a function of the wall inertia parameter; spatial profiles of real and imaginary parts of the eigenfunctions at neutral stability for (c)mode-(i)(ˆρ=0.631); (d) mode-(ii) ( ˆρ=12.73); (e) mode-(iii) ( ˆρ=21.02). Here,ˆ T0=5, ˆpext =2.98 and Re =50. and 10 as a function of the wall inertia parameter ˆρ; at neutral stability the eigenfunction profile of the oscillatory mode has two (three) extrema in the real (complex) part of thewallprofile(figure 15c), meaning the number of extrema can change over a period of oscillation for these walls of finite thickness. For this choice of parameters the primary oscillatory mode of the lower branch (mode-(i)) is stable for ˆρ=0, becoming unstable as the wall inertia parameter, ˆρ, increases (figure 15a), while the corresponding oscillation frequency decreases (figure 15b). The perturbation growth rate for lower branch mode-(i) exhibits a local maximum at ˆρ≈10 before asymptoting toward zero as the wall inertia parameter continues to increase. Hence, this mode of instability approaches stability with decreasing oscillation frequency as the wall gets heavier. However, as the wall inertia parameter increases a second mode of oscillation also becomes unstable at ˆρ≈12.72 (figure 15a) with larger frequency than mode-(i) (figure 15b); we term this mode-(ii), which also has a perturbation wall profile with two extrema (figure 15d) albeit with a narrow boundary layer at the upstream end of the profile. Unlike the primary mode, the growth rate of this instability continues to increase as ˆρincreases for these parameter values, while the corresponding oscillation frequency again approaches zero (figure 15b). As the wall mass parameter becomes even larger, we eventually observe another mode becoming destabilised for ˆρ≈21.02 with larger frequency (figure 15b) which we term mode-(iii), where the wall profile again exhibits two extrema with a narrow upstream boundary layer (see eigenfunction wall profile in figure 15e). Further increases in the wall inertia parameter destabilises mode-(iv) (figure 15a,b). Note that, in accordance with our naming convention, the oscillation 934 A28-25 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press M.A. Herrada, S. Blanco-Trejo, J. Eggers and P.S. Stewart PONCE-TORRES,A.,RUBIO,M.,HERRADA, M.A., EGGERS,J.&MONTANERO, J.M. 2020 Influence of the surface viscous stress on the pinch-off of free surfaces loaded with nearly-inviscid surfactants. Sci. Rep. 10 (1), 16065. RYZHAKOV, P.B., ROSSI,R.,IDELSOHN,S.R.&ONATE, E. 2020 A monolithic Lagrangian approach for fluid–structure interaction problems. Comput. Mech. 46, 883–899. SHAPIRO, A.H. 1977 Steady flow in collapsible tubes. Trans. ASME J. Biomech. Engng 99 (3), 126–147. SHTERN, V. 2018 Models of fold-related hysteresis. Phys. Fluids 30 (5), 2–7. SHTERN,V.&HUSSAIN, F. 1996 Hysteresis in swirling jets. J. Fluid Mech. 309, 1–44. SNOEIJER, J.H., PANDEY,A.,HERRADA,M.A.&EGGERS, J. 2020 The relationship between viscoelasticity and elasticity. Proc. R. Soc. Lond. A476, 20200419. STEPHANOFF,K.,PEDLEY, T.J., LAWRENCE,C.&SECOMB, T.W. 1983 Fluid flow along a channel with an asymmetric oscillating constriction. Nature 305, 692–695. STEWART, P.S. 2010 Flows in flexible channels and airways. PhD thesis, University of Nottingham. STEWART, P.S. 2017 Instabilities in flexible channel flow with large external pressure. J. Fluid Mech. 825, 922–960. STEWART, P.S., HEIL,M.,WATERS, S.L. & JENSEN, O.E. 2010 Sloshing and slamming oscillations in a collapsible channel flow. J. Fluid Mech. 662, 288–319. STEWART, P.S., WATERS, S.L. & JENSEN, O.E. 2009 Local and global instabilities of flow in a flexible-walled channel. Eur. J. Mech. (B/Fluids) 28 (4), 541–557. WALTERS, M.C., HEIL,M.&WHITTAKER, R.J. 2018 The effect of wall inertia on high-frequency instabilities of flow through an elastic-walled tube. Q. J. Mech. Appl. Maths 71 (1), 47–77. WANG,D.,LUO,X.Y.&STEWART, P.S. 2021aEnergy analysis of collapsible channel flow with a nonlinear fluid-beam model. J. Fluid Mech. 926,A2. WANG,D.,LUO,X.Y.&STEWART, P.S. 2021bMultiple steady and oscillatory solutions in a collapsible channel flow. Intl J. Appl. Mech. 13 (4), 2150058. WHITTAKER, R.J., HEIL,M.,JENSEN, O.E. & WATERS, S.L. 2010 Predicting the onset of high-frequency self-excited oscillations in elastic-walled tubes. Proc. R. Soc. Lond. A466 (2124), 3635–3657. XU,F.,BILLINGHAM,J.&JENSEN, O.E. 2013 Divergence-driven oscillations in a flexible-channel flow with fixed upstream flux. J. Fluid Mech. 723, 706–733. XU,F.,BILLINGHAM,J.&JENSEN, O.E. 2014 Resonance-driven oscillations in a flexible-channel flow with fixed upstream flux and a long downstream rigid segment. J. Fluid Mech. 746, 368–404. XU,F.&JENSEN, O.E. 2015 A low-order model for slamming in a flexible-channel flow. Q. J. Mech. Appl. Maths 68 (3), 299–319. ZHANG,S.,LUO,X.Y.&CAI, Z. 2018 Three-dimensional flows in a hyperelastic vessel under external pressure. Biomech. Model. Mechanobiol. 17 (4), 1187–1207. 934 A28-32 https://doi.org/10.1017/jfm.2021.1131 Published online by Cambridge University Press