Effect of base-state curvature on self-excited high-frequency oscillations in flow through an elastic-walled channel

We investigate the development of flow-induced small-amplitude high-frequency oscillations occurring in a rigid fluid-conveying channel with a section of the upper wall replaced by a taut elastic sheet. We consider an axially nonuniform base state, caused by a negative transmural pressure (internal minus external), which results in the sheet bending inward, and then examine the evolution of oscillatory perturbations about this state. Due to the curvature of the base state, any normal displacement of the sheet will lead to a change in the axial stretching of the sheet, which is an order of magnitude higher than would be the case for perturbations to a uniform sheet. This stretching provides an additional restoring force that can be dominant. We derive a modified tube law to describe the wall mechanics with this additional effect and combine it with an existing fluid model to obtain a complete description of the system. At leading order, we obtain a one-dimensional eigenvalue problem for the frequencies and mode shapes of the oscillatory perturbations. The nonlinear interaction between the perturbations and the base-state curvature manifests itself as an additional integral term in the eigenvalue problem, corresponding to the total axial stretching in the wall. The normal modes are neutrally stable at leading order, and their slow growth or decay is determined by considering the global energy budget of the system. The stability of the system can be expressed in terms of a critical Reynolds number of the mean flow through the channel. We explain the behavior of the system as two key dimensionless parameters F and K are varied. These quantify the dimensionless axial tension and base-state curvature effects, respectively. Numerical simulations are used to identify distinct flow regimes in parameter space, which we explore further in detail using asymptotic analysis. This reveals that at large K, i.e., strong base-state curvature effects, the leading-order axial profiles of the modes are forced to adjust in order to eliminate the need for significant stretching in the wall. This has the effect of stabilizing both the fundamental and, to a lesser extent, the first harmonic modes, which results in the first harmonic becoming the most unstable mode.

Experimental studies of such systems are usually conducted using a Starling resistor setup. This consists of a length of elastic tubing connected between two rigid pipes placed within a pressure chamber. A flow is driven though the pipes and tubing by either applying a pressure difference between the two ends or using a volumetric pump at one end. When the transmural pressure (internal minus external pressure) across the tube wall becomes sufficiently negative the wall buckles and becomes susceptible to large changes in shape in response to small changes in pressure. A closely related physical model consists of a rigid two-dimensional (2D) channel with a segment of one wall replaced by an elastic membrane under axial tension [8] (see Fig. 1). This system has attracted considerable theoretical attention due to its relative simplicity in comparison with the Starling resistors while still experiencing phenomena such as self-excited oscillations.
Experiments demonstrate the strong interaction between the elastic tube and the fluid it contains. A wide variety of behaviors have been observed ranging from high-frequency "flutter" to lowfrequency "milking" [9]. The different types of oscillation form in response to different boundary conditions, mean flows, external pressures, tube-wall properties, and fluid properties [10,11]. Much experimental work has focused on the onset and saturation of self-excited large-amplitude oscillations in the Starling resistor. In this case, the cross-sectional area of the tube remains near its maximal value for the majority of the period and oscillations are shown to develop when the flow rate exceeds a critical value [10][11][12].
Simplified theoretical analyses of such flows have been based on spatially one-dimensional models [13] and require an ad hoc description of flow dissipation and tube-wall mechanics [14]. Despite their relative simplicity, such models do incorporate the basic elements required to describe self-excited oscillations and are able to capture many of the qualitative behaviors observed in experiments. Numerical simulations of the two-dimensional channel analogy reveal a wide variety of behavior [15]. This has assisted in the development of more rational theoretical models.
When the axial tension in the elastic wall is relatively low, the transmural pressure causes the (unstable) steady solution to have a highly deformed wall shape. Theoretical analysis is complicated by the nonuniform base state and complex flow structures in the large-amplitude oscillations, however the underlying mechanisms behind such instabilities are now beginning to be fully understood through recent work by Stewart et al. [14,16] and Xu et al. [17,18].
Alternatively, when the wall stiffness and/or tension is larger, instabilities can occur about less deformed static states, with higher-frequency oscillations. Jensen and Heil [13] introduced a model that could describe such an instability in a 2D channel in which one wall is replaced by an elastic membrane. In this parameter regime, the oscillatory fluid flow induced by the wall motion typically has a two-region structure comprising an inviscid flow in the core, bounded by thin viscous Stokes layers near the channel walls. Conservation of volume implies that, during phases when the wall moves inward, some of the fluid in the collapsible section must be displaced into the upstream and 113901-2 downstream rigid channels. The oscillatory motion of the wall therefore generates axial sloshing flows that increase (decrease) the volume flux in the downstream rigid section while decreasing (increasing) the volume flux in the upstream rigid section when the wall moves inward (outward).
The energy fluxes into and out of the system are the kinetic energy flux through the tube ends, the work done by pressure and viscous forces at the tube ends and at the wall, and the losses due to dissipation. If an instability is to grow, the total energy of the system, on average, needs to increase over time.
The mechanism of sloshing instabilities identified by [13] for 2D channel flow relies on an increase in kinetic energy into the system and involves interactions between the high-frequency axial sloshing flow and the boundary conditions at the tube ends. This provides a means for energy transfer from the mean flow to a growing oscillatory perturbation. The boundary conditions at the upstream and downstream ends determine how this axial oscillatory flow is partitioned into upstream and downstream components. Any oscillatory perturbation to the steady flow at either end will increase the kinetic energy flux there. If the upstream and downstream conditions are such that the oscillatory perturbation to the axial flow is greater at the upstream end than the downstream end, then there will be a net increase of kinetic energy into the system.
The essential elements of this sloshing instability mechanism are also present in threedimensional tube flow [19,20]. In a sufficiently long tube, the cross-sectional flows are in general relatively unimportant and the axial sloshing flows can still dominate the dynamics. The sloshing mechanism relies on changes in the cross-sectional area. The efficient extraction of energy from the mean flow is then promoted when tube oscillations are performed about a nonaxisymmetric mean state. (Small-amplitude deformations of a circular cross section result in area changes that are quadratic in displacement amplitude. The induced sloshing flows are then an order of magnitude smaller, resulting in a smaller influx of kinetic energy and a weaker instability.) Hence instabilities are more likely to occur in tubes with noncircular cross sections. Such arrangements are also more representative of the buckled configuration about which oscillations typically develop in experiments [21]. However, it is noted that the high-tension sloshing mechanism requires largeramplitude oscillations at the upstream end of the flexible tube. Most experimental studies to date have observed larger oscillations at the downstream end. The sloshing mechanism is therefore probably not relevant to the lower-tension regimes that have been explored experimentally, but this does not preclude it being important at higher axial tensions.
Previous work on high-frequency oscillations in a 3D tube, occurring about a relatively undeformed base state, focused on the effect of axial tension relative to azimuthal bending in the tube wall [22]. As a function of the axial coordinate z ∈ [0, 1], the amplitude of the nth axial mode of the transverse wall perturbations was found to resemble the graph of sin[(n + 1/2)π z], with the n = 1 (axial-mode-1) perturbation always being the most unstable. Intuitively it might be thought that the sloshing instability described above would always operate most efficiently with an axial-mode-1 oscillation since all of the displaced fluid is forced out of the ends of the tube after each contraction. For n > 1, some fluid will oscillate entirely within the elastic-walled section of the tube rather than helping to force an oscillatory flow at the ends, thus making no contribution to the kinetic energy flux at the ends. However, if the tube wall oscillates about a more deformed initial configuration, the axial-mode-1 oscillations require significantly more axial stretching relative to the mode-2 oscillations. In the latter case, axial displacements of the wall can eliminate the need for such stretching (see Fig. 2). Since stretching is much harder than bending in thin walls, the restoring forces in a mode involving significant stretching will be expected to be larger than for those that do not. We therefore anticipate that either the axial-mode-1 modes oscillate at a higher frequency in response to these larger forces or the axial profile of the mode alters to compensate for the nonuniform base state. For sufficiently buckled and thin elastic walls, either change may provide sufficient stabilization to the fundamental mode for the first harmonic to become the most unstable.
In this paper we examine the effect of nonuniform initial wall configurations on the onset of high-frequency small-amplitude self-excited oscillations in a rectangular channel. This work represents a proof-of-concept theoretical consideration of this additional physical effect. The simplified geometry and linearized model allow analytical progress to be made while still retaining the essential physics. The aim is to provide insights into and understanding of the effect of base-state curvature on the stability of the system, which can then inform the development of more realistic models.
We extend the theoretical work of Whittaker et al. [22] to investigate the interaction between the curvature of the initial configuration and axial tension in the wall, and the effects this has on the stability of perturbations of the wall. In Sec. II we define the problem and derive a model that describes the motion of the elastic wall in response to pressure in the fluid. We combine this with the description of the fluid presented by Whittaker et al. [19,20] to obtain a model that determines the profile and frequency of perturbations. We then consider the energy budget of the system to obtain a stability criterion that determines the slow growth or decay of the oscillating disturbances over time. In Sec. III we outline the numerical methods used to solve the model system. The results of these calculations, evaluated in different regions of parameter space, are presented in Sec. IV and compared to asymptotic approximations and reduced models derived in Appendixes A-C. The key results are summarized and discussed in Sec. V. Concluding remarks can be found in Sec. VI.

II. MODEL
We consider the flow of a fluid through a rigid, rectangular channel of width b, height a, and length L, aligned with the coordinate axes (x * , y * , z * ), respectively, as shown in Fig. 3 of the upper wall at y * = a, between z * = z * 1 and z * = z * 2 , is replaced by an elastic sheet with density σ and wall thickness h that is allowed to move freely. At z * = z * 1 and z * = z * 2 the elastic section of the wall is clamped to the rigid sections fixing the position and axial slope of the elastic sheet.
We assume that the sheet, in its undeformed state, is subject to a uniform in-plane stress of magnitude T 0 /b in the axial direction, where T 0 is the extensional force applied at the ends of the undeformed sheet. The material of the sheet is assumed to be linearly elastic for small-amplitude deformations, with Young's modulus E and Poisson's ratio ν. We define the bending stiffness K and extensional stiffness D as We also assume that the elastic wall is connected to external springs, each with spring constant κ, that are uniformly distributed with number density n per unit area and aligned perpendicular to the surface of the undeformed sheet [see Fig. 3(a)]. Choosing nκ = O(K/a 4 ) allows us to consider this setup as a two-dimensional analog of flow though an elliptical or partially collapsed circular cylindrical elastic tube [8], where the restoring force created by the springs is analogous to the restoring force from azimuthal bending in the three-dimensional case. (As shown in [23], in a 3D system with an initially elliptical cross section, a single azimuthal mode dominates deformations of the wall, resulting in a restoring force per unit area from azimuthal bending that scales as K/a 3 , i.e., K/a 4 times the magnitude a of the deformation.) The fluid inside the channel has density ρ and viscosity μ. We impose a steady axial volume flux Q at the downstream end of the channel and a fixed fluid pressure p * = p * up at the upstream end. We apply no-slip boundary conditions on the upper (moving) and lower channel walls and periodic boundary conditions on the sidewalls. An external pressure p * ext acts on the top of the sheet. Deformations to the undeformed state are then caused by the transmural pressure p * − p * ext . We consider the case where p * ext is steady and causes an initial steady deformation of maximum amplitude d = δa to the elastic portion of the wall. We term this steady deformed configuration the base state and then consider oscillatory perturbations of typical amplitude e = a about this state. The amplitudes e and can be time dependent if the oscillations are growing or decaying.
The timescale of the oscillations turns out to be which arises from a balance between oscillatory components of the axial fluid inertia (with the inertial pressure scale ρL 2 /τ 2 ; see Sec. II C) and the force from the compression and extension of the springs (a surrogate for azimuthal bending, with pressure scale anκ ∼ K/a 3 ). We assume that any growth or decay of the oscillations, i.e., variation in , takes place over a timescale that is much longer than τ . Thus time derivatives of may be neglected relative to those of the oscillatory functions. The formal asymptotics of this approach were considered in [24], where it was found to be appropriate for this system.

A. Dimensionless groups and parameter regime
There are a number of dimensionless groups in the problem, the most important of which are as follows. First we have several geometric parameters which describe the aspect ratio, positions of the flexible wall ends, wall thickness, base-state deformation amplitude, and oscillatory perturbation amplitude, respectively.

113901-5
Then there are three parameters related to the wall mechanics: The parameter F gives the ratio of restoring forces from axial tension to those from the springs, while M gives the ratio of wall inertia to fluid inertia; K turns out to be the appropriate measure of the strength of the nonlinear effect from the initial wall curvature. [The factor of (z * 2 − z * 1 ) is included in K to simplify later calculations.] Finally, there are two parameters related to the fluid mechanics: where τ is the timescale defined in (2) above. The Womersley number α gives the relative importance of unsteady inertia to viscous effects. The Strouhal number St gives the relative importance of unsteady to convective inertia. For this initial study, we consider a regime in which This means we have a long thin-walled tube with negligible wall inertia, which has a slightly deformed base state. The tube wall then undergoes high-frequency oscillations with an amplitude much smaller than that of the base-state deformation. The small amplitude of the base state is a largely a mathematical convenience, but is acceptable for a proof-of-concept exploration of the effects of base-state curvature. At least near the onset of instability, the oscillations will be small in amplitude, so δ is justified here.

B. Model for the wall mechanics
Changes in the transmural pressure result in the wall being displaced from its undeformed position. The deformation is resisted by the wall's elasticity, by the wall's inertia, and by the force from the springs. In a two-dimensional deformation, the wall elasticity acts through axial bending and axial stretching. In the long-wavelength limit considered here, the resistance to axial bending is small compared with the other forces [22]. Therefore, we will neglect axial bending effects. Wall inertia is included in the derivation below, but will later be neglected, following our assumption M 1 above. We describe the deformation by its components η * = η * (z * ) in the vertical direction and ζ * = ζ * (z * ) in the axial direction [see Fig. 3(b)]. Under the deformation, a short element of wall, of undeformed length dz * , is stretched to a length The wall displacements are written in terms of their base-state and oscillatory perturbation components (denoted by overbars and tildes, respectively): Recall that the parameters δ 1 give the relative sizes of the base state and oscillatory perturbations, soη andη are expected to be O (1). Moreover, we must have max |η| = 1 so that the deformation amplitude in the base state is precisely d = δa. The scales in the expressions for ζ * follow from geometrical considerations and to obtain a leading-order balance between ζ * z * and η * 2 z * in (9). The nondimensional cross-sectional area of the channel is the O(δ 3 ) correction term in the cross-sectional area arises due to the horizontal displacement of the wall. The transmural pressure, written aš is decomposed into steady and oscillatory components and nondimensionalized using the scale K/a 3 associated with bending of the elastic walľ The steady external pressure is then defined as The leading-order tension in (9) is nondimensionalized using (17) and (18) to obtain For later use, energy and energy fluxes are nondimensionalized on respectively, based on the kinetic energy and kinetic energy fluxes in the steady flow. Now, using (17)-(23), Eq. (13) for the normal displacements is nondimensionalized, yielding, when = 0 and to O( ), respectively, and where k 0 = a 4 nκ/K = O(1) and k 2 = a/b = O(1) represent the strength of the springs relative to azimuthal bending and the aspect ratio of the channel; k 0 and k 2 are analogous to the O(1) constants employed by Whittaker et al. [22] for flow through an elastic cylindrical tube that depend on the geometry of its undeformed state. Similarly, nondimensionalizing Eq. (16) for the axial displacements, we find, when = 0 and to O( ), respectively, and where the dimensionless parameters M and K are as defined in Sec. II A.

Fixing the axial tension in the base state
Deforming the elastic wall into its base-state configuration while keeping the end points fixed would increase the axial tension in the wall. To concentrate on the effect that base-state curvature has on the dynamics of the wall, we allow a small movement of an end point in order to leave the axial tension in the deformed base state unchanged from the initial value of T 0 /b. Any increase (decrease) in the axial tension is then as a result of the lengthening (contracting) of the wall as it oscillates. This allows comparison of perturbations about different base states to be considered independently of changes in the axial tension. To achieve this, we assume that when the wall is deformed into its base-state configuration, the edge at (y, z) = (1, z 2 ) is moved horizontally closer to (y, z) = (1, z 1 ) so that With a small-amplitude base-state deformation (δ 1), Eq. (27) can be integrated to givē for some constant c. Then, from (23), we must have c = 0 to ensure that (29) holds. Hencē Integrating (31) and applyingζ = 0 at z = z 1 , we obtain which sets the required horizontal displacement at (y, z) = (1, z 2 ).

Tube law for the wall displacements
The oscillatory components of the horizontal and vertical displacements are fixed at the boundaries, i.e.,ζ = 0 andη = 0 at z = z 1 and z = z 2 . Integrating (28) assuming M 1 and δ 1, we obtainζ where the constant of integration has been found by integrating again by parts and applying the boundary conditions onζ andη given above. Using (31) and (33) to eliminateζ z andζ z from (25) and (26) and assuming M 1, we obtain These equations govern how the base stateη and oscillationsη are forced by the transmural pressuresP andP. The perturbations are subject to the boundary conditionsη =η = 0 at z = z 1 and z = z 2 . Recall also the normalization max |η| = 1 of the base-state deformation.
Equation (34a) for the initial steady deformation takes the form of a standard tube law, with the transmural pressureP balanced by the restoring forces k 2 Fη zz from the axial tension and k 0η from the springs. Equation (34b) for the subsequent perturbations also contains an additional term arising due to the nonlinear interaction between the perturbations and the base-state curvature. Physically, the integral corresponds to the total axial stretching (increase in length) of the wall, and the overall term is the product of resulting axial tension perturbations and the base-state curvature.

C. Model for the fluid
The fluid velocity u * and pressure p * within the channel are described by the Navier-Stokes equations for an incompressible Newtonian fluid: Integrating the continuity equation ∇ * · u * = 0 over a cross section of the channel, we have whereẑ is the unit vector in the axial (z * ) direction and A(z * ) is the cross-sectional area of the channel. A mean flow is driven along the channel by imposing a steady axial volume flux at the downstream end and fixing the pressure at the upstream end, i.e., where Q is the prescribed fluid flux. [Other boundary conditions may be used instead, e.g., a flux condition or a pressure condition at the upstream and downstream ends of the channel. However, in the analysis below we focus on the boundary conditions imposed by (37) as these ensure that the instability mechanism described in the Introduction is at its most powerful. Prescribing the flow rate at the downstream end ensures that no kinetic energy is lost there.]

Nondimensionalization
Following Whittaker et al. [19], the fluid velocity is decomposed as u * ≡ u * ⊥ + w * ẑ into axial and transverse components. Then u * ⊥ , w * , and p * are each nondimensionalized and decomposed into a steady, oscillatory, and nonlinear component. Here we focus on the steady and oscillatory components (denoted by overbars and tildes) and write Due to the slender geometry of the channel, the axial and transverse components of the velocity are nondimensionalized on different scales. The leading-order axial velocityw is nondimensionalized using the velocity scale U = Q/ab of the background flow, while the leading-order transverse oscillatory flow is nondimensionalized using the velocity scale of the wall oscillations. The scales for other two velocity components follow from the continuity equation. The steady and oscillatory components of the pressure are nondimensionalized by viscous and inertial scales, respectively, relative to the appropriate velocities.
Equations (35)-(37) are nondimensionalized using (17) and (38). The O(1) components of (35) and (37) govern the steady component of the fluid flow. These are simply the steady Navier-Stokes equations, written in a dimensionless form. The solution of a steady high-Reynolds-number flow problem is not straightforward; however, the exact details of the steady flow are unnecessary for evaluating the leading-order oscillatory flow and associated energy transfer. All the information that is required is contained in the (nondimensional) time-averaged volume flux [19], which is Q from (37b). The oscillatory component of the flow is determined by considering (35) at O( ). In particular, the axial component of this yields The dynamics of the fluid are parameterized by the Womersley number α, the Strouhal number St, and the aspect ratio , as defined in Sec. II A. The (steady) Reynolds number is given by For our model, we need a description of the fluid flow in response to periodic oscillations of the elastic wall, valid in the long-wavelength ( 1) high-frequency (α, St 1) small-amplitude (δ, 1) regime specified in (6). In this parameter regime the O( ) axial component of the the Navier-Stokes equation (39) reduces tow Because of the high frequency of the oscillations (α 1), the oscillatory axial velocity has a plug flow profile in the core, with passive Stokes layers adjacent to the elastic wall. The magnitude of the plug flow is determined by conservation of mass, given by (36). An axial pressure gradient is required to accelerate the oscillating fluid. The transverse component of (35) gives ∇ ⊥p = 0 to O( ), implying that transverse components of the pressure and axial velocity are uniform in the cross section in this parameter regime. Combining (19) and (36), to O( ) we find We seek solutions for the oscillatory plug flow component of the form where ω is oscillation frequency, and the circumflexed variables (ŵ,p,η) are the (possibly complex) amplitudes of the leading-order oscillations as functions of axial position. Equations (41) and (42) then become which then combine to givep providing a relationship between oscillatory axial pressure and the vertical displacement of the wall.

113901-11
Substituting the ansatz (43) into the boundary conditions (37) and making use of (44), we find the following boundary conditions on the oscillatory amplitudep of pressure perturbation: Physically, the fixed pressure at the upstream end (z = 0) means there must be no oscillatory pressure perturbation there. The fixed flux at the downstream end (z = 1) means that there must be no oscillatory flux perturbation there, which translates to no oscillatory component to the axial pressure gradient.

D. Combining the solid and fluid models via the transmural pressure
To combine the solid mechanics model of Sec. II B and the fluid mechanics model of Sec. II C we need to relate the dimensionless transmural pressuresP andP to the dimensionless fluid pressures p andp. Starting from P * = p * − p * ext and substituting from (21), (22), and (38c), we have Equating the leading-order steady and oscillatory terms and using (5) and (2), we find

Base-state deformation
To simplify the calculations, we now make the additional assumption δα 2 St 1. From (48), this then means that the steady transmural pressureP is dominated by the applied external pressurē p ext , and we can neglect the contribution from the steady viscous pressure dropp. Hence we havē (This assumption is just a mathematical convenience to allow solutions forη to be written down explicitly. It should not have any significant qualitative effect on the results.) We now take the steady transmural pressureP to be negative and uniform in z and then solve (34a) for the base-state deformationη, subject toη(z 1 ) =η(z 1 ) = 0 and min{η} = −1. We find where γ 2 = k 0 /k 2 F, with the dimensionless steady transmural pressure taking the valuē Using (49) and (22) to relate toP to p * ext , we then have that the dimensionless amplitude of the base-state deformation is given by When the axial tension is small relative to the spring restoring force, i.e., F 1, or equivalently γ 1, the base-state displacement profile and amplitude are approximated bȳ Conversely, when the axial tension is large relative to the restoring force from the springs, i.e., F 1 or γ 1, the base-state displacement profile and amplitude are approximated bȳ

Oscillatory perturbation
Combining (34b), (45), and (48), we obtain a single equation forp(z), valid for z 1 < z < z 2 , where is a measure of the total nonlinear stretching, and hence the increase in axial tension, due to the oscillatory perturbation in the elastic portion of the wall. As we will see below, the size of I relative to the other terms in (57a) and relative top zz in (57b) will be important in determining the character of the solutions to the system. In the rigid sections of the wall (0 < z < z 1 and z 2 < z < 1) (57) is replaced with the constraint η = 0; combining with (45), this then impliesp We apply the boundary conditions (46) at z = 0, 1 to obtain Matching between the different regions at z = z 1 and z 2 requires the continuity of pressure and axial volume flux, i.e.,p andp z are continuous. Equation (57) is then subject to the boundary conditions The system of (57) and (60) defines an eigenvalue problem forp(z) and ω given a base stateη(z) and parameters F and K. We expect to find a countably infinite set of normal modes. These modes are purely oscillatory at leading order, with any growth or decay being a first-order effect.

E. Stability criterion and growth rate
Following Whittaker et al. [22], we now determine the slow growth or decay of the normal modes by considering the energy budget of the system. Multiplying (35) by u * , integrating over the volume of the channel, and averaging over one period of wall oscillations, the energy budget of the fluid, nondimensionalized on the scales (24), is partitioned as [22] where E is the mean rate of working by the fluid on the wall, T is the mean flux of kinetic energy through the tube ends, S is the mean rate of working by pressure forces at the tube ends due to the oscillatory flow, and D is the mean rate of dissipation by the oscillatory flow. The final term in (61) accounts for the amplitude of oscillations varying slowly over time and involves the period-averaged  kinetic energy E f of the oscillatory flow. The energy fluxes may then be decomposed into two components: the energy flux arising from the mean flow alone and energy flux arising through interactions of the mean flow with the oscillatory flow. Averaging over one oscillation of the wall, the former are shown to cancel out from (61). The latter contribution is expressed as [22] T = 3 The mean kinetic energy of the oscillatory flow is dominated by the axial plug flowŵ and is given by We now consider the energy budget of the elastic wall, integrated over the whole wall and averaged over one period of oscillation. The only source of energy is the work done E by the fluid on the wall and the only store of energy is elastic energy E s . (There is no kinetic energy since inertia is neglected in the wall, M 1.) Using the scalings (24), we therefore have where the flux E is obtained from the energy budget of the fluid (61) and the mean elastic energy is [22] wherep † denotes the complex conjugate. Eliminating E between (61) and (65), we find Writing (62), (63), and (65) in terms ofp using (44a) and (45), we find that and The latter implies that the time-averaged energy of the mode is partitioned equally between the kinetic energy in the fluid and the elastic energy of the wall. Both (67) and (68) are proportional to 2 . We therefore have which suggests that the amplitude of the oscillations grows or decays exponentially, with the growth rate being given by allowing (70) to be written in the form = 0 (Re − Re c ), where

III. NUMERICAL METHOD
The eigenvalue problem defined by (57) and (60) was solved numerically using a spectral method. The functionp(z) was approximated as a sum of N Chebyshev polynomials and the governing equations were evaluated at the Chebyshev collocation points k j = cos[(2 j − 1)π/2N] for j = 1, . . . , N, yielding an eigenvalue problem of the form for some N × N matrices A and B. This was then solved using the eig in-built MATLAB eigenvalue problem solver. The integral I in (57) was obtained by evaluating the sum of the product ofη zzpzz and the Clenshaw-Curtis weights at each of the N collocation points [25]. The spectral accuracy of the approximation was verified by comparison with solutions obtained when k 2 = K = 0. The accuracy of solutions is not bounded by the numerical scheme but by the number N of Chebyshev polynomials discretizing the domain; increasing N reduces the error in the solution. When N > 11 the error in the frequency remains smaller than 10 −8 . Solutions to (57) and (60) with K = 0 were also compared with previously published results by Whittaker et al. [22]. Once solutions forp had been obtained, the growth rates and critical Reynolds numbers Re c were evaluated using (70) and (71).

IV. RESULTS
The numerical simulations presented below were conducted with z 1 = 0.1, z 2 = 0.9, k 0 = 11.074 87, and k 2 = 1.704 41 for consistency with results presented by Whittaker et al. [22]. Unless stated otherwise, solutions are normalized such that the pressure perturbationp satisfies The dynamics of the wall is governed by the eigenvalue problem defined by (57) and (60) and is dependent on the shape of the base state (50) (examples of which are illustrated in Fig. 4 for different F). The system is parametrized by two key parameters: F, the strength of axial tension relative to azimuthal bending, and K, a measure of the nonlinear stretching effects as a result of the base-state curvature in the wall. For each point in (F, K) parameter space we find a set of normalmode solutions and oscillation frequencies. We define the primary and secondary modes as the solutions to (57) and (60) resulting in the smallest and subsequent smallest oscillation frequencies ω, respectively, for a given base state and parameter values. Previous work in high-frequency oscillations about an undeformed base state found that the primary and secondary modes adopt axial-mode-1 and axial-mode-2 shapes, respectively, and that the primary mode is always the most unstable. Our key result is that in some (large K) regions of parameter space, each mode increases its axial wave number to add an extra turning point. As a result of this, the secondary mode can become the most unstable. This is discussed in more detail in the following sections.  Figure 5(a) illustrates the oscillation frequency of the primary and secondary modes for varying F and K. We observe that increasing K results in an increased frequency of oscillations. The frequency of the primary mode increases the most, but it remains strictly below that of the secondary mode. Figure 5(a) also allows us to identify three distinct parameter regimes: F 1, F = O(1), and F 1. When F 1, ω appears to be independent of K to leading order. As F becomes O(1) we observe that increasing K results in faster oscillations. When F is large, collapse of the scaled oscillations [illustrated in Fig. 5(b)] suggests that the dynamics of the wall become dependent on a single parameter F = F/K. Figure 6(a) shows the corresponding values of Re c illustrating the Reynolds number above which the wall oscillations grow exponentially. We observe that for K above some critical value, there is a finite window in F where Re c is smaller for the secondary mode than the primary mode. This means that within this window the system can be unstable to the secondary mode while the primary-mode oscillations decay over time. Figure 6(b) illustrates this region of (F, K) parameter space. We see that for sufficiently large F there exists a critical value of K above which the secondary-mode perturbations lose stability before the primary mode. We define K c (F ) as the value of K at which Re c corresponding to the primary and secondary modes are equal; this is minimized by F ≈ 0.29 where K c ≈ 1.50. Figure 6(b) also illustrates that in the large F limit the stability curve asymptotes to F c ≈ 0.37. Figure 7 shows how the total nonlinear axial stretching I varies as a function of F and K for the primary mode. The magnitude of I will be important in determining the character of the solutions in different regions of parameter space. We see that I → 0 as F → 0, and there is collapse for I as a function of K/F when this ratio is large.
We now look at the pressure profiles and the shape of the elastic wall in each of the parameter regimes identified above. We begin by addressing the limit of F 1 before considering F = O(1) and F ∼ K 1. Simulations will be compared to analytic approximations and to solutions of reduced models (illustrated by dashed or dotted curves in Figs. [8][9][10][11] in Appendixes A-C.

A. Small axial tension: F 1 and K F 3/2
Small F represents the case where the axial tension in the undeformed wall is weak relative to the spring restoring force. Figure 4 shows that for O(γ −1 ), where γ 2 = k 0 /Fk 2 . The oscillatory componentη of the wall displacement shares a similar three-region structure [see Fig. 8(b)]; however, this feature is not observed inp [ Fig. 8(a)].
In line with the standard nomenclature within the theory of matched asymptotic expansions, we will refer to the two small regions near z = z 1 and z = z 2 as boundary layers. This terminology should not be confused with any "boundary layers" (such as viscous Stokes layers or Womersley layers) that may occur in the fluid close to the channel walls.
The dynamics of the elastic wall are dominated by the interior of the domain where there is little curvatureη zz (Fig. 4). As a result, the nonlinear axial stretching I [as defined in (57b)] is dominated by contributions from the boundary layers near z = z 1 and z = z 2 . The term in I is found to make a negligible contribution to (57a), which explains the lack of K dependence inη in the bulk region connecting the boundary layers (see Fig. 8) and in the oscillation frequency ω of the wall for sufficiently small F [ Fig. 5(a)]. Figure 7 illustrates how the magnitude of I varies with respect to F and K. The collapse of the curves in Fig. 7(b) suggests that for sufficiently large γ 2 K, we have I = O((γ 2 K) −1 ). This observation will be exploited in the analysis below. for reference, to illustrate the effect curvature in the base state has on the structure ofη. We see that in the boundary layers, the curvature in the base state results in a change in sign of the oscillatory component of the wall displacement:η decreases to a local minimum before increasing with a large gradient. As F increases, the length scale of the boundary layers becomes comparable to that of the interior of the domain, the linear scaling of I with γ −2 breaks down, and we enter the F = O(1) regime.
Motivated by the numerical results, we use a matched asymptotic expansion to further explore this parameter regime. Recalling γ 2 = k 0 /k 2 F, we introduce a scaled oscillation frequencyω 2 = ω 2 /k 0 and a scaled nonlinear stretching λ = (Kγ 2 /k 0 )I = O(1). The governing equation (57) then becomes where We consider expansions in inverse powers of γ : In Appendix A we find a matched asymptotic solution (A34) that satisfies (74) to O(γ −1 ) and gives values for the composite functions and constants that form the solution. The solution is obtained by considering O(γ −1 ) boundary layers near z = z 1 and z = z 2 that are connected by an O(1) bulk interior region away from the boundaries. We obtain a composite solution of the form wherep II 0 (z) andp II 1 (z) are, respectively, the leading-order and first-order correction pressure profiles obtained from the bulk interior region and are given by (A18).
We see that the leading-order pressure profile in (77) is determined by the dynamics of the wall in the bulk of the domain and is independent of the axial tension F, and the base-state curvature effects K, which appear only in the terms γ −2 = Fk 2 /k 0 and λ = (K/k 0 )γ 2 I = O(1). Decreasing γ (increasing F) then has little effect on the amplitude of the oscillations in the bulk region, but instead shortens the width of the bulk region. The base-state curvature effects contribute to the pressure profile at O(γ −1 ) and have an effect in the boundary layers close to z = z 1 and z = z 2 only [entering (77) in the second and third lines].
The composite solution (77) satisfies (74), but is dependent on the constants λ i , i 0. To determine these constants (77) must be combined with (76b) and (75). The leading-order term on the left-hand side of (75) must balance with the highest-order nonzero term on the right-hand side. To obtain the numerically observed scaling I = O((Kγ 2 ) −1 ), we find that (77) is constrained by the conditions that K F 3/2 and that the leading-order normal-mode solution in the two boundary layers must make no combined contribution to the integrals on the right-hand side of (75). This implies that when K F 3/2 , the base-state curvature effects in the boundary layers are strong enough to forceη to be orthogonal toη zz there. This makes (77) independent of K to leading order. Full details can be found in Appendix A.
The leading-order solution is the solution to the full eigenvalue problem with K = 0. Corrections are found sequentially by computing I from the solution at the previous order with (57b) and then using it as a forcing term in (57a). The full details can be found in Appendix B. When F = O(1) and K 1, the base-state curvature effects dominate over the restoring force from the springs. For sufficiently large K, Eqs. (57) and (60) are satisfied by solutions of the form The leading-order solution is the solution of an eigenvalue problem in which KI appears as an unknown O(1) constant. Since this means I must be smaller than an order-of magnitude estimate using (57b) would suggest, we must have i.e., the solution forη 0 = −p 0zz /ω is forced to be orthogonal to the base-state curvature. The constant KI is chosen to achieve this. At the next order, there is a similar eigenvalue problem with another constant that must be set to satisfy an integral constraint on the solution. The full details can be found in Appendix B.
The leading-order terms in (78) and (79) describe the dynamics of the wall when K is sufficiently small or large, respectively, and provide useful bounds on the oscillation frequency ω and critical Reynolds number Re c for intermediate K [see Figs. 9(a) and 9(b)]. In both limits, the leading-order solutions are independent of K. Figure 10 illustrates the pressure-perturbation and wall-displacement profiles for the primary and secondary modes, for the case F = 1, with K = 10 −0.5 and K = 10 2.5 (red and blue curves, respectively). These parameter values are marked by crosses in Fig. 9 and show typical behaviors in the K 1 and K 1 regimes. As K increases from the K 1 regime to the K 1 regime the local maximum inη shifts downstream. Additionally, as a consequence ofη being forced to be orthogonal toη zz at leading order, each mode gains an additional half oscillation in z. The primary mode therefore resembles an axial-mode-2 mode, the secondary resembles an axial-mode-3 mode, etc.
The sloshing mechanism responsible for driving instabilities to grow works most efficiently with odd axial modes rather than even axial modes, since less of the displaced fluid is trapped between pairs of opposing extrema. This has the effect, as K increases, of increasing the stability of the primary mode resulting in the secondary mode becoming the most unstable. Figure 9(b) shows that the switching of the order in which the primary and secondary modes become unstable occurs in the intermediate regime where K = O(1). However, considering the transition of the system from the K 1 to the K 1 asymptotics provides insight into the dynamics of the system at this critical point.

C. Large axial tension: F 1 and F /K O(1)
When F 1 the axial tension in the undeformed wall is large relative to the spring restoring force. The base state is approximated by (55) and has uniform curvature along the axis of the channel. Figure 5(b) suggests that for sufficiently large F and K it may be possible to describe the system in terms of a single parameter F = F/K, with the oscillation frequency ω scaling as K 1/2 .
We therefore write ω = ω/K 1/2 = O(1) and find in Appendix C that the dynamics of the system can be approximated by (C4) and (C3b). These equations define a leading-order eigenvalue problem for ω = ω (F ), valid when F 1 and F/K O(1). The system is solved numerically, and some example solutions for the pressure and wall displacement amplitudes can be found in Fig. 11. The region F 1 and F/K O(1) includes part of the boundary in (F, K) space on which the order in which primary-and secondary-mode solutions lose stability switches. As K varies, we see a transition similar to that shown in Fig. 10 (and discussed above in Sec. IV B) for F = O (1).
The solutions to the reduced system (C4) and (C3b) for F 1 allow us to identify the asymptotic critical value F c of F = F/K at which the order of stability loss switches between the primary and secondary modes. We find that F c ≈ 0.37. This is shown by the dashed line in Fig. 6(b). The critical solutions for the primary and secondary modes that have the same Re c are illustrated by the red curves in Fig. 11. At this critical point, the wall displacements corresponding to the primary and secondary modes resemble the graphs of − sin(5π z/2 − π/4) and sin(5π z/2 − π/4), respectively.

V. DISCUSSION
We explain the behavior of the system as the two key dimensionless parameters K and F are varied. A schematic of the parameter space is shown in Fig. 12. The parameter F measures the strength of the axial tension in the wall relative to the restoring force from the springs (which is a surrogate for the azimuthal bending effects that would be present in a 3D tube). By considering the origins of the terms in (57), the parameter K can be interpreted as measuring the strength of the basestate curvature effects relative to the restoring force from the springs. When K is sufficiently small, the base-state curvature effects are then a small perturbation to the K = 0 solution. For sufficiently large K, the base-state curvature effects dominate. These two regimes are separated by a transition region occupying K = O(F 3/2 ) when F 1 and K = O(F ) when F 1. Inspecting (57), we would anticipate that for large enough K, either the K term dominates and the eigenfrequencies ω become O(K) or we have I = O(K −1 ) to prevent the K term from becoming too large. Our numerical results show that the latter option occurs. This may be explained on theoretical grounds (at least for lower-axial-order modes) as follows. With a balance between the K and ω 2 terms in (57),p would be forced to be proportional to the base-state curvatureη zz . Boundary layers would then be required to satisfy the boundary conditions onp. However, in all parameter regimes, the possible boundary-layer solutions lack sufficient modes that decay into the interior and would result in the forcing of short-wavelength oscillatory behavior throughout the domain. Thus loweraxial-order modes cannot have I = O(1) for K F. To consider the behavior in more detail it is useful to split the cases F O(1) and F 1, as the latter will involve boundary layers.
First we will consider F O (1). In this case, the dominant balance in (57) when K = 0 involves terms of O(F ). So for K F the K term is a small perturbation to that system. At leading order, we still have the K = 0 solution, which can be used to compute I from (57b). This in turn can be reinserted into (57a) to induce the first-order correction.
As K increases to O(F ), the size of the perturbation grows to become as large as the leading-order solution, and we enter the transition region. For K F, we must have I = O(F/K) 1, so from (57b) the shape of the wall oscillations (η, which is proportional top zz ) must be almost orthogonal to the base-state curvatureη zz . Writing the solution as an asymptotic expansion in K, the leading-order solutionp 0 must satisfy z 2 z 1η zzp0zz dz = 0. Then I must be chosen in (57a) to ensure that this is the case. In effect, the value of KI in (57a) is a Lagrange multiplier for the constraint of (almost) no axial stretching.
As K increases from O(F ) and leaves the transition region, each of the modes must adjust to become more orthogonal to the base-state curvature. This results in each mode gaining an additional half-wave oscillation in z. So the fundamental (primary) mode becomes axial-mode 2, the first harmonic (secondary mode) becomes axial-mode 3, etc. This has the effect of reducing the instability of the fundamental mode relative to the first harmonic, as axial-mode-2 modes give less sloshing flow at the upstream end, as more displaced fluid can be transferred between the two opposing extrema. When K exceeds a critical value K c (F ) the most unstable mode switches from the fundamental to the first harmonic. The region where the stabilities are switched is sketched in Fig. 12 and shown quantitatively in Fig. 6(b).
When F 1, the situation is complicated due to the presence of boundary layers of width O(F 1/2 ) in both the base state and the oscillatory perturbations. In this case the base-state curvature is negligible outside the boundary layers, and hence the base-state curvature effects do not directly affect the bulk oscillatory solution in the interior (although there is some influence through the matching conditions from the boundary layers). Effects similar to those described for F O(1) occur in the boundary layers, with the solution being a perturbation to the K = 0 solution for K F 3/2 , and I appearing in (57a) as a Lagrange multiplier for K F 3/2 . However, because there is no significant K effect on the mode shapes, the stability switching does not occur for F 1.

113901-24
EFFECT OF BASE-STATE CURVATURE ON SELF-…

VI. CONCLUSION AND FUTURE WORK
This study of flow through a channel with one surface replaced by an elastic sheet has demonstrated the importance of base-state curvature on the dynamics of the wall in subsequent self-excited oscillations. Building on the results of Whittaker et al. [22], which considered perturbations about a uniform initial state in a 3D tube, we derived a system that governs the base-state and oscillatory components of the wall displacement valid for long-wavelength high-frequency small-amplitude displacements. To simplify matters in this initial study, we considered flow in a 2D channel, as a surrogate for the full 3D system.
The base-state deformation is modeled by (50) when a uniform external pressure is applied. Subsequent small-amplitude oscillations about this configuration are governed by the eigenvalue problem (57) and (60), which is solved numerically to find a countable set of normal modes with different frequencies. Considering the energy budget of the system predicts that these modes grow or decay with a rate given by (70). The stability criterion (71), expressed in the form of a critical mean-flow Reynolds number, determines when each mode is stable.
The model system depends on two key nondimensional parameters: F, a measure of the axial tension, and K, a measure of the base-state deformation. The main finding of this work is that for sufficiently large F and K (see Fig. 12), the effect of the base-state curvature causes the mode with the second-lowest frequency to be the most unstable, ahead of the fundamental mode.
The model has limitations, for example, the 2D setup, neglecting the effects of wall inertia relative to inertia in the fluid, neglecting bending in the elastic wall relative to tension effects, and neglecting the contribution to the transmural pressure from the steady viscous pressure drop along the channel. However, we anticipate that the key physical processes described here will be applicable to flows in more realistic practical applications. Care also needs to be taken when applying these results, noting that they are only applicable to an elastic wall undergoing small-amplitude long-wavelength high-frequency oscillations.
In this paper we restricted our attention to a channel with fixed cross-sectional aspect ratio k 2 = 1.704 41 and fixed strength of springs relative to azimuthal bending k 0 = 11.074 87, but chose to survey the full (F, K) parameter space. Typically, experimental realizations of self-excited oscillations utilize a Starling resistor, consisting of an axisymmetric elastic tube, that buckles nonaxisymmetrically under a sufficiently negative transmural pressure. To fully understand the mechanisms of such phenomena, the analysis presented may need be extended to consider the full 3D system. This includes consideration of full Shell equations for the wall and nonlinear effects in the fluid. There may also be a need to improve the fluid model to incorporate the effects of flow past a nonuniform base state. A further area of interest would be to investigate oscillations about a larger-amplitude initial base-state configuration, i.e., δ = O(1). In the derivation of the wall equations (25)-(28), terms that are quadratic and higher order in the horizontal and vertical deformations will be significant at leading order. This will lead to increased complexity in the governing equations; however, it may provide important insight into the dynamics of the system at later times when the amplitude of instabilities has grown to become comparable to the radius of the tube.

ACKNOWLEDGMENT
The authors would like to acknowledge the financial support from the Engineering and Physical Sciences Research Council Grant No. EP/N007212/1 to undertake this work. When the axial tension in the elastic sheet is small (F 1), the base stateη is approximated by (53), an example of which is illustrated in Fig. 4. The system (57) is then approximated by Below we derive a matched asymptotic solution to (A1), subject to the boundary conditions (60), for γ 1, i.e., F 1. This matches solutions obtained in the boundary layers near z = z 1 (region I) and z = z 2 (region III) over a length O(γ −1 ) [see Figs. 8(c) and 8(d)], together with the bulk interior region (region II) of length O(1) [see Fig. 8(b)]. (Superscripts denote solutions in each of the three regions.) The only contributions to I come from the two boundary layers. It is convenient to write I = I I + I III , where I I and I III are the contributions from regions I and III, respectively.

Boundary layer near z = z 1 (region I)
In the boundary layer close to z = z 1 we define a scaled coordinatez = γ (z − z 1 ). We then setp =p I (z) and seek solutions that are normalized such thatp Iz (0) = γ −1 , i.e.,p z (z 1 ) = 1. Equation (A1) then becomesp with and is subject top To match with region II we must also suppress the exponentially growing solutions to the homogeneous equation. Together with (A5) this gives four boundary conditions on the solution at each order. We seek solutions of the form Substituting these into (A3), it follows that to leading order and O(γ −1 ), respectively, Enforcing (A5) and ensuring thatp I 0 andp I 1 do not grow exponentially large asz → ∞, we find To O(γ −2 ) and O(γ −3 ), respectively, Eq. (A3) becomes which, with (A5) and suppressing exponentially growing terms, are satisfied bỹ The region I solution (A6a) may be substituted into (A4) to obtain an expansion for I I of the form 113901-26 where Writing (A6a) in terms of the O(1) variable z, the pressure in the boundary layer near z = z 1 is approximated bŷ We note that we retain the exponentially decaying terms in (A13), despite them being O(γ −2 ), since they contribute to the leading-order wall displacement, evaluated using (45). The unknown constantsω 0 ,ω 1 , λ 0 , and λ 1 will be determined later through the matching of (A13) with solutions obtained in the bulk of the domain (region II) and the outer boundary layer (region III).

Bulk region (region II)
In the region connecting the boundary layers near z = z 1 and z = z 2 , we have z = O(1) and p =p II (z) and the dominant balance in (A1) iŝ p II zz +ω 2pII = O(γ −2 ).
For the numerically observed scaling I = O((γ 2 K) −1 ) to hold we need λ 0 = O(1). The leadingorder term on the left-hand side of (A38) must balance with the largest nonzero term on the righthand side. Hence we must have γ 3 K O(1). If γ 3 K 1, i.e., K F 3/2 , then to balance λ 0 with an O(1) term, we must equate it with a higher-order term on the right-hand side of (A38) and require I 0 = 0. The larger γ 3 K is, the further along the series on the right-hand side of (A38) one must go to match the λ 0 term.
In the region of parameter space where γ 3 K 1, we must have I 0 = I I 0 + I III 0 = 0. Writing this in terms of the normal-mode solution, via (45) and (57b), we thus have Consequently, when γ 3 K 1, base-state curvature effects in the boundary layers are strong enough to force the leading-order normal-mode displacements to be orthogonal to the base-state curvature, hence making (A34) independent of K to leading order.

APPENDIX B: REDUCED MODELS FOR MODERATE AXIAL TENSION F = O(1)
When the axial tension F is O(1), increasing K (the base-state curvature effect) results in a faster frequency of wall oscillations. Figure 9 suggests the existence of K 1 and K 1 regimes, in both of which the dynamics of the wall is described by a system independent of K. In Appendix B 1 and B 2, we formulate reduced models that capture the dynamics of wall when K 1 and K 1, respectively. subject top 1zz = 0, z 1p 1z −p 0 = 0 at z = z 1 , (B11a) p 1zz = 0,p 1z = 0 at z = z 2 .