Multiple Steady and Oscillatory Solutions in a Collapsible Channel Flow

We study flow driven through a finite-length planar rigid channel by a fixed upstream flux, where a segment of one wall is replaced by a pre-stressed elastic beam subject to uniform external pressure. The steady and unsteady systems are solved using a finite element method. Previous studies have shown that the system can exhibit three steady states for some parameters (termed the upper, intermediate and lower steady branches, respectively). Of these, the intermediate branch is always unstable while the upper and lower steady branches can (independently) become unstable to self-excited oscillations. We show that for some parameter combinations the system is unstable to both upper and lower branch oscillations simultaneously. However, we show that these two instabilities eventually merge together for large enough Reynolds numbers, exhibiting a nonlinear limit cycle which retains characteristics of both the upper and lower branches of oscillations. Furthermore, we show that increasing the beam pre-tension suppresses the region of multiple steady states but preserves the onset of oscillations. Conversely, increasing the beam thickness (a proxy for increasing bending stiffness) suppresses both multiple steady states and the onset of oscillations.


Introduction
Fluid flow through flexible conduits in the human body can exhibit a wide variety of interesting physiological phenomena [Heil and Hazel, 2011]; such flows can be investigated experimentally using a Starling Resistor, where fluid is driven through a segment of externally pressurised flexible tubing by either a prescribed driving pressure or a prescribed upstream flux [e.g. Bertram, 1982Bertram, , 1986. In particular, pressure driven flow can exhibit a phenomenon known as 'flow limitation', where the flow rate along the tube does not continue to increase as the driving increases: the large flow speeds collapse the tube through the Bernoulli effect inhibiting the flow [Bertram and Castles, 1999;Bertram and Tscherry, 2006]. Similarly, flux driven systems can exhibit the associated phenomenon of 'pressure drop limitation', where the pressure difference between the ends of the collapsible tube does not increase with increases in flow rate [Bertram, 1986;Bertram et al., 1990]. In physiology, flow limitation can occur during tidal breathing, when the flow rate of air being expelled from the lung becomes maximised during forced expiration [Tantucci, 2013].
In some cases these highly collapsed steady configurations can co-exist with other, more inflated, configurations of the vessel. Such multiplicity in steady solutions is evident in Starling Resistor experiments, where the system can exhibit hysteresis i.e. different steady configurations at the same operating point dependent on the history and termed an 'open-to-closed transition' [Bertram et al., 1991;Bertram and Castles, 1999]. Such co-existing steady states are also evident in theoretical models of flow in collapsible tubes [Hazel and Heil, 2003;Heil and Boyle, 2010]. In particular, lumped parameter models indicate that the collapsible tube can typically exhibit three co-existing branches of steady solutions across a range of parameters [Armitstead et al., 1996]; these three branches have been termed the upper branch (where the tube is mostly inflated), the lower branch (where the tube is highly collapsed) with an intermediate branch between them. The upper (lower) steady solutions is connected to the intermediate branch through a saddle-node bifurcation at the upper (lower) branch limit point; the upper and lower steady branches are stable to non-oscillatory perturbations and the intermediate branch is always unstable [Armitstead et al., 1996].
In physiology, flow limitation in the lung airways is often accompanied by wheezing, associated with rapid fluttering of the airway wall [Grotberg and Gavriely, 1989]. Similarly, Starling Resistor experiments investigating flow limitation can sometimes exhibit spontaneous self-excited oscillations which fall into distinct frequency bands [Bertram and Pedley, 1982;Bertram et al., 1990Bertram et al., , 1991. Such selfexcited oscillations are also evident in theoretical models of flow in collapsible tubes, such as lumped parameter models [Bertram and Pedley, 1982;Armitstead et al., 1996], cross-sectionally averaged one-dimensional models [Jensen, 1990] and full three-dimensional models [Heil and Boyle, 2010;Whittaker et al., 2010].
In this study we consider a theoretical model for a planar analog of the Starling Resistor experiment, formed by removing a segment of one wall of a rigid channel and inserting an elastic wall. This planar analog exhibits multiple steady solutions in some parameter regimes [Luo and Pedley, 2000;Heil, 2004;Stewart, 2017;Wang et al., 2021;Herrada et al., 2021], with a three branch structure qualitatively similar to the collapsible tube models. In addition, this collapsible channel system exhibits transition to self-excited oscillations from both the lower (collapsed) branch of steady solutions [Heil, 2004;Stewart, 2017] as well as the upper (inflated) branch of steady solutions [Herrada et al., 2021;Wang et al., 2021]. Fully developed upper branch oscillations exhibit an upstream propagating hump along the compliant segment; this wave is suppressed by reflection at the upper rigid segment and replaced by a new wave at the downstream end of the compliant segment [Wang et al., 2021]. By contrast, fully developed lower branch oscillations exhibit a constriction at the downstream end of the compliant segment which propagates up and downstream over a narrow range [Luo and Pedley, 1996;Wang et al., 2021].
In particular, we use a model where the wall takes the form of a thin prestressed (massless) elastic beam with resistance to both bending and stretching (Sec. 2) [Cai and Luo, 2003;Luo et al., 2008;Wang et al., 2021]. We use this model to explore how the size of the region of multiple steady states is influenced by the system parameters, particularly the beam pre-tension and bending stiffness (Sec. 3). Finally, we also use this model to explore the onset of self-excited oscillations from both the upper and lower branches of steady solutions, showing that in some cases these two families of oscillations merge together (Sec. 4).

The Model
We consider an incompressible Newtonian fluid (of density ρ and viscosity µ) driven through a finite-length rigid channel with uniform width D. We consider a parabolic inlet flow driven through the channel with flux Q against a fixed outlet pressure p 0 . One segment of the upper wall is replaced by a plane strained elastic beam of thickness h subject to a uniform external pressure p e . The corresponding lengths of the upstream and downstream rigid and elastic segments are denoted L u , L d and L, respectively. We parameterize the domain using a Cartesian coordinate system with origin at the intersection between the upstream rigid segment and the compliant segment on the entirely rigid wall (see Fig. 1). Time is denoted t.
We model the flexible wall as a massless elastic beam and denote the axial pre-tension along the beam as T . The extensional and bending stiffnesses of the beam are denoted as EA and EJ, respectively. Here E is the Young's modulus of the material while A and J are the cross-sectional area and the second moment of inertia of the cross-sectional area of the beam, respectively. The undeformed elastic beam is parameterized by the coordinate l (where 0 l L).

Governing Equations
To form non-dimensional variables we scale all lengths on D, velocities on Q/D, time on Q/D 2 and pressures on the inertial scale according to where p (p) denotes the dimensional (dimensionless) fluid pressure. The corresponding dimensionless parameters take the form wherec λ ,c κ andT are the dimensionless extensional, bending stiffnesses and pretension of the elastic beam, respectively,Re is the Reynolds number andp e is the external pressure. The corresponding dimensionless lengths of the upstream, downstream and collapsible segments of the channel and the dimensionless beam thickness are scaled as respectively. We henceforth focus on the dimensionless quantities and drop the tildes for simplicity.

Fluid Equations
The governing equations for the (two-dimensional) fluid are the incompressible Navier-Stokes equations in the form, where u = (u, v) is the planar fluid velocity and the Newtonian stress tensor σ takes the form where I is the identity matrix and the superscript T represents the matrix transpose.

Beam Equations
We denote the two components of beam deformation as x b = (x b (l, t), y b (l, t)) in terms of the beam coordinate l; using a modified constitutive law for the massless beam, the dimensionless governing equations for the beam can be written as (details see Wang et al. [2021]) where θ is the angle between the rigid wall and the tangent vector of the deformed beam (see Fig. 1), κ is the curvature of the beam and λ is the principal stretch of the beam, which can each be computed in terms of the beam deformation as In addition, σ 1 and σ 2 are the tangent and normal components of the fluid traction on the beam, wheret andn are the tangent and normal unit vectors of the deformed beam (see Fig. 1).

Boundary Conditions
We prescribe a parabolic inlet flow with unit flux in the form We assume the no-slip condition along the rigid walls as well as continuity of velocity between the elastic beam and the fluid in the form The two ends of the elastic beam are fixed to the rigid wall in the form

The Finite Element Method
A finite-element method is used to solve the coupled fluid-beam system. We divide the fluid domain into three sections, denoted as A, B and C for the upstream, compliant and downstream compartments ( Fig. 1), respectively [Luo and Pedley, 1996;Luo et al., 2008]. In section B, we use an adaptive mesh that consists of rotational spines that connect fixed nodes in the bottom wall with nodes in the deformable beam [Cai and Luo, 2003]. Then nodes are seeded along these spines covering region B. Each spine can rotate around its fixed node on the rigid wall. Hence, nodes in section B can move along the rotational spine as the beam is deformed. Further details of the numerical method are provide elsewhere [Luo et al., 2008;Hao et al., 2016]. A mesh of 36657 elements is used for the numerical solutions in this study with time-step ∆t = 0.01. Convergence tests of grid-and time-stepindependence were carried out between three meshes and two choices of time-step (for details see Wang et al. [2021]).

Parameter Choices
In our previous studies we modeled the elasticity of the beam by varying the extensional stiffness c λ [Luo et al., 2008;Wang et al., 2021], which is proportional to the bending stiffness of the beam since c κ = h 2 /12 c λ , and fixed the other beam parameters including the pre-tension (T = 0), beam thickness (h = 0.01) and external pressure (p e = 1.95). Conversely, in this study we fix c λ = 1600 and examine the role of these other parameters on both the steady (Sec. 3) and unsteady (Sec. 4) behaviour of the system. In particular we note that increasing the beam thickness allows us to alter the bending stiffness of the structure independently of the extensional stiffness. Following Luo et al. [2008], we fix the geometry of the channel according to L u = L = 5, L d = 30.

Multiple Steady Solutions
We discuss the predictions of the steady system considering the role of increasing pre-tension (Sec. 3.1) and beam thickness (Sec. 3.2).

Role of Pre-tension
In order to investigate the role of pre-tension T in setting the steady beam shape, Fig. 2 summarises the steady solutions of the maximal and minimal beam positions ( Fig. 2a), several typical beam shapes (Fig. 2b) and the upper and lower limit points ( Fig. 2c) with different values of pre-tension with fixed wall thickness h = 0.01. In particular, the maximal (y max , solid line) and minimal (y min , dash-dot line) beam deflection on the y−direction as a function of Reynolds number Re is shown in Fig.  2(a) for T = 0 (black), T = 1 (red) and T = 5 (blue) at fixed external pressure p e = 1.1. For low Reynolds numbers the beam is fully inflated (i.e. y min = 1); this is the upper branch of steady solutions. As the Reynolds number increases the elastic beam becomes increasingly collapsed though the Bernoulli effect. Notably, for T = 0 the system enters a region with three steady states for 354.56 Re 362.78 similar to our previous model [Wang et al., 2021]. As the Reynolds number further increases, the system again exhibits a unique steady state for Re > 362.78; this branch is the lower branch of steady solutions. The upper and lower branches of steady solutions are connected by an intermediate branch, which they intersect at the upper and lower limit points, respectively. A similar region with multiple steady states is observed for T = 1 (355.79 Re 358.75), although the beam pre-tension has significantly narrowed the range of Reynolds numbers for which it is evident for this choice of external pressure. The region with multiple steady states vanishes entirely for large pre-tension (T = 5, blue lines in Fig. 2a).
To illustrate these steady configurations in detail, Fig. 2 vertical line in Fig. 2a), the steady wall shape is fully inflated and is unique (solid line in Fig. 2b). Conversely, at Re = 360, located in the region with multiple steady states for T = 0 (see the dashed vertical line in Fig. 2a), there are three possible steady wall shapes (Fig. 2b). Whereas for T = 5 (Fig. 2c), the system has a unique steady state for both Reynolds numbers and the wall is significantly less deflected due to the larger pre-tension. Fig. 2(d-f) summarises the region with multiple steady states for T = 0 ( Fig.  2d), T = 0.1 (Fig. 2e) and T = 1 (Fig. 2f) in parameter space spanned by the external pressure and Reynolds number. In each case we see a triangular region with multiple steady states with external pressure below a critical value and Reynolds number above a critical value. As T increases, this critical point is gradually pushed to larger Reynolds number and lower external pressures. However, the region with multiple steady states (Fig. 2f)   the system by shifting the critical point for multiple solutions across the parameter space.

Role of Beam Thickness
To investigate the role of increased bending stiffness on the steady configuration of the beam, Fig. 3 demonstrates the steady solutions with various beam thickness, plotting the maximal (y max ) and minimal (y min ) beam deflection against the Reynolds number Re for three beam thicknesses (h = 0.01, 0.086, and 0.122) with fixed pre-tension T = 0 and external pressure p e = 1.1 in Fig. 3(a). Similar to Fig.  2, the system shows multiple steady states for h = 0.01 (354.56 Re 362.78); this region of multiple solutions narrows with increased beam thickness h = 0.086 (354.87 Re 355.8). Additionally, the steady system is unique for sufficiently large beam thickness h = 0.122. Fig. 3(b-d) summarises the corresponding upper and lower limit points in parameter space spanned by the external pressure p e and Reynolds number Re for the same three beam thicknesses, showing that the region with multiple steady states narrows as h increases, but the critical point for the existence of multiple solutions does not vary much across the parameters tested. In summary, this figure demonstrates that increasing the beam thickness (and hence the bending stiffness) eventually suppresses multiple steady states; however, this is achieved by narrowing the tongue while holding the critical point fixed.

Self-excited Oscillations
In order to test the stability of the system for a given parameter combination we apply a small increment to the steady solution (here we use the steady solution corresponding to a 1% increase in the extensional stiffness c λ ). The system is deemed stable if the unsteady solution converges to its corresponding steady solution, and unstable if the perturbation grows [Drazin, 2002]. The state between stable and unstable is termed as neutrally stable. In this model the perturbation grows in an oscillatory manner and saturates into a finite amplitude limit cycle. In this study, we report the dynamics of this oscillatory limit cycle and ignore the initial transient.
In particular, we present phase portraits of the oscillation as a function of the wall pressure at the upstream and downstream ends of the compliant segment (e.g. Fig.  5c-f, i-l and Fig. 6b-e) and compute the fluid pressure on the wall at the channel midpoint time-averaged over a period of oscillation (e.g. Fig. 5a,g, Fig. 6a, Fig. 8 and Fig. 9). In this section, we focus on the unsteady solutions for fixed extensional stiffness (c λ = 1600), testing the stability of both the upper and lower steady branches. In particular, we investigate oscillatory solutions with no pre-tension and low beam thickness (Sec. 4.1), large beam thickness (Sec. 4.2) and large pre-tension (Sec. 4.3).

Multiple Oscillatory Solutions with Low Bending Stiffness
We first analyse unsteady simulations of perturbations to the upper and lower steady branches with no pre-tension (T = 0) and a thin beam (h = 0.01). Here we extend our previous analysis [Wang et al., 2021] to examine the parameter space spanned by the external pressure and Reynolds numbers, as shown in Fig. 4   Wang et al.
[2021] we held the external pressure p e = 1.95 throughout). The steady behaviour of the system is similar to the cases presented in Fig. 2(d-f) and 3(bd), where three steady states can co-exist in a narrow tongue between the upper and lower limit points (blue shaded region marked by filled symbols). Testing the stability of these steady solutions we find that both the upper and lower steady branches can independently become unstable to oscillations (each via a supercritical Hopf bifurcation) in the neighbourhood of the region of parameter space with multiple steady states. The computed neutrally stable points from the upper and lower steady branches are marked as filled circles, and the neutral stability curves are denoted as dot-dashed lines connecting these neutral points. In this parameter space the upper branch is unstable within a narrow tongue to the left of the trace of the upper branch limit points. This upper branch neutral curve intersects the trace of the upper branch limit points at a co-dimension 2 (Takens-Bogdanov) point [Glendinning, 1994;Strogatz, 2018] at p e ≈ 2.25, Re ≈ 174.5 (we term this the upper Takens-Bogdanov point). For external pressures less than this upper Takens-Bogdanov point the upper branch oscillation is eventually restabilised as the Reynolds number increases via an interaction between the oscillatory limit cycle and the upper branch limit point [see details in Wang et al., 2021]. However, this unstable tongue extends to slightly larger external pressures than those which admit see multiple steady states. For external pressures greater than the upper Takens-Bogdanov point the upper branch oscillation instead restabilises via a second Hopf bifurcation as the Reynolds number increases. This interaction is explored further in Fig. 5, where we examine the time-averaged midpoint pressure of fully developed upper branch oscillations (Fig. 5a), their corresponding period (Fig. 5b) and illustrate several limit cycles of oscillation ( Fig. 5c-f). As in Wang et al. [2021], this upper branch oscillation is associated with an overall decrease in the time-averaged midpoint pressure compared to the steady state (Fig. 5a), and the period of oscillation increases as the upper limit point is approached (Fig. 5b). The oscillation develops a rather complicated limit cycle (Fig. 5c-f), similar to the upper branch oscillations reported in Wang et al. [2021].
Furthermore, the lower branch of steady solutions is unstable within a tongue to the right of the region with multiple steady states. This tongue is stabilised for large external pressure (p e 2.12), and tracks close to the curve of lower branch limit points, merging at a second (lower) Takens-Bogdanov point at Re ≈ 266.8, p e ≈ 1.48. The neighbourhood of this lower Takens-Bogdanov point is explored in Fig. 5(g-l), plotting the time-averaged midpoint pressure (Fig. 5g), period of oscillation (Fig. 5h) and several limit cycles (Fig. 5i-l) for fixed external pressure p e = 1.48. During oscillation the midpoint pressure is again less than the corresponding steady configuration. Just beyond the critical Reynolds number for lower branch instability, the oscillation grows and saturates into an elliptical limit cycle of larger period than the upper branch oscillations (open circles in Fig. 5h) but of much smaller amplitude (open circles in Fig. 5g). A phase portrait for point L1 is shown in Fig.  5(i). For Reynolds numbers sufficiently close to critical (Re 266.8) this smallamplitude limit cycle is maintained over the lifetime of our numerical simulations. However, for larger Reynolds numbers (Re 270.3) the system visits this limit cycle only transiently and eventually saturates into a limit cycle of larger amplitude (open triangles in Fig. 5g) and shorter period (open triangles in Fig. 5h). These nonlinear oscillations are a direct continuation of the oscillations bifurcating from the upper branch steady state, displaying analogous phase portraits (Fig. 5j,k). Furthermore, across the range of Reynolds numbers which exhibit the small-amplitude limit cycle growing from the lower steady branch (266.8 Re 270), it turns out that the system also exhibits a saturated limit cycle analogous to the upper branch oscillations for these points as well (amplitude and period shown as open triangles in Fig. 5(g,h) and a typical phase portrait is shown in Fig. 5(l). Hence, this system appears to exhibit hysteresis in the oscillatory behaviour across a small range of parameters with two possible branches of fully developed oscillations (open circles and open triangles in Fig. 5h); these two branches merge together at Re ≈ 270.3 and then continue to large Reynolds numbers (open triangles in Fig. 5h). There are two possibilities to explain this observation: either the system exhibits two coexisting limit cycles across a range of parameters, or these reported lower branch (small-amplitude) limit cycles are in fact long transients which eventually grow and saturate to a limit cycle analogous to the upper branch oscillations (note that we see no evidence of this transition for simulations close to the critical Reynolds number over the long time intervals considered in our simulations). Either way, for large enough Reynolds numbers the system exhibits a merged family of oscillations (open triangles in Fig. 5h) which grow from the lower steady branch but are a direct continuation of nonlinear oscillations bifurcating from the upper steady branch.
It emerges that for external pressures less than the lower Takens-Bogdanov point, this merged family of oscillations becomes more prominent with a typical example for p e = 1.1 shown in Fig. 6, showing the bifurcation diagram (Fig. 6a) and four typical limit cycles in the neighbourhood of the region with multiple steady states ( Fig. 6b-e). As the Reynolds number increases the system becomes unstable along the upper branch of steady solutions, exhibiting a decrease in the time-averaged midpoint pressure with an elaborate limit cycle around the upper branch steady state (Fig. 6b). As the Reynolds number increases into the region with multiple steady states the oscillation encompass all three steady solutions (Fig. 6c, d), and there is now no interaction with the intermediate steady state (a prominent feature of purely upper branch oscillations, as discussed in Wang et al. [2021]). As the Reynolds number continues to increase into the region with only a lower steady state, the oscillation exhibits a qualitatively similar limit cycle (Fig. 6e).
To further explore the nature of these merged oscillations, Fig. 7 examines the oscillation in more detail at a point in parameter space which exhibits multiple steady states, showing time-trace of the wall pressure at the two ends of the compliant segment (Fig. 7a) and several snapshots of the flow-field and pressure contours in the channel (Fig. 7b-e). The fully developed (nonlinear) oscillation retains many of the characteristics of the upper and lower branch oscillations from which it has merged (these individual oscillations are discussed in Wang et al. [2021]). Similar to the upper branch oscillations, the limit cycle involves an upstream propagating hump in the wall profile which is suppressed by interaction with the upstream rigid segment and replaced by a new upstream propagating hump originating at the downstream end of the compliant segment. However, like the lower branch oscillations, the wall profile exhibits significant indentation into the channel (and so much lower fluid pressures). This greater wall indentation leads to the formation of a low pressure region near the end of the compliant segment (Fig. 7d), which creates a (weak) vorticity wave propagating along the downstream rigid segment (Fig. 7e, f), similar to oscillations about a collapsed steady state [Luo and Pedley, 1996]. Inter- estingly, the wall pressures at the two ends of the compliant segment are much closer to being in phase than the upper or lower branch oscillations alone. In summary, the figure illustrates the flow profile of our new family of self-excited oscillations, showing that it retains characteristics of both the upper and lower oscillations that have been reported elsewhere [Wang et al., 2021].

Oscillations with Large Beam Thickness
To evaluate the role of the beam thickness, h, in the unsteady response of the system, Fig. 8 presents bifurcation diagrams showing the onset of self-excited oscillations, plotting the time-averaged midpoint pressure alongside the corresponding steady midpoint pressure as a function of the Reynolds number for three beam thicknesses (h = 0.086, 0.094 and 0.012) with zero pre-tension and fixed external pressure (p e = 1.1). For beam thickness h = 0.086 (Fig. 8a), the system exhibits multiple steady states in the range 354.87 Re 355.8. In this case the oscillations from the upper and lower steady branches merge into one family of oscillations across the region with multiple steady states (similar to Fig. 6). When the beam thickness is increased to h = 0.094, the multiple steady states vanish but the system still exhibits transition to self-excited oscillations at Re ≈ 344.07 (denoted as the dashed line), which grow in amplitude as the Reynolds number increases (Fig. 8b). However, as the beam thickness is further increased to h = 0.122 the system no longer becomes unstable to oscillations and the steady system has a unique solution for all external pressures (Fig. 8b). In summary, increasing the thickness of the beam (and hence the bending stiffness) suppresses the onset of self-excited oscillations as well as the region with multiple steady states.

Oscillations with Large Pre-tension
In order to evaluate the role of pre-tension T in the onset of self-excited oscillations, Fig. 9 Fig. 9. Bifurcation diagrams of midpoint wall pressure p mid as a function of the Reynolds number Re at h = 0.01, pe = 1.1 for (a) T = 1 and (b) T = 5. The steady midpoint pressure p mid is denoted as solid (stable) and dashed (unstable) lines. The time-averaged midpoint wall pressure p avg mid is denoted as filled triangles, the upper and lower branch limit points are denoted as squares.
compared to the steady midpoint pressure as a function of the Reynolds number for large pre-tension T = 1 (Fig. 9a) and T = 5 (Fig. 9b) for fixed external pressure (p e = 1.1) and beam thickness (h = 0.01). For T = 1, the steady system exhibits multiple steady solutions (Fig. 2) and the oscillations initiated close to the upper and lower steady branches merge into one family of oscillations across the region with multiple steady solutions (Fig. 9a). However, for much larger pre-tension (T = 5) the system exhibits a unique steady state for all Reynolds numbers but this still becomes unstable to self-excited oscillations at Re ≈ 338.3. As before, the time-averaged midpoint pressure is lower than the corresponding steady value and decreases with the Reynolds number (Fig. 9b). Therefore, the onset of oscillations is preserved with increasing pre-tension despite the loss of multiple steady states.

Discussion
In this study we revisit a theoretical model for flow in a planar collapsible channel where the flexible wall is modelled as a pre-stressed elastic beam using a modified nonlinear constitutive law [Wang et al., 2021], investigating the influence of the external pressure, beam pre-tension and thickness (a proxy for bending stiffness) on the steady and unsteady behaviour of the system. The model was solved numerically using the finite element method. Similar to the Starling Resistor experiments [Bertram et al., 1990;Bertram and Castles, 1999] and previous models of flow in collapsible channels and tubes [Armitstead et al., 1996;Heil, 2004;Stewart, 2017], our model predicts that the steady system can exhibit multiple co-existing states, consisting of upper, intermediate and lower steady branches. The model predicts that this region of multiple steady states is suppressed by increasing either the wall pre-stress (Fig. 2) or the wall bending stiffness (Fig. 3). In the former the critical point for multiple steady states is postponed to larger Reynolds numbers and lower external pressures as the pre-stress increases (Fig. 2d-f), whereas in the latter the critical point does not move significantly but the region with multiple steady states narrows as the bending stiffness increases (Fig. 3b-d).
Previous studies have indicated that self-excited oscillations can (independently) grow from either the upper and the lower branches of steady solutions [Wang et al., 2021] in the neighbourhood of the region with multiple steady states (Fig. 4, 5). Our model predicts that these two families of oscillations eventually merge into a single family of oscillatory behaviour for sufficiently low external pressure (Fig. 6). This new merged family of oscillations retains characteristics of both the upper and lower branch oscillations (Fig. 7). This new family of oscillations is preserved as the beam pre-stress increases, despite the loss of multiple steady states (Fig. 9). Conversely, this merged family of oscillations is suppressed by increasing the beam thickness (i.e. increasing the bending stiffness of the beam, Fig. 8).