Linear and weakly nonlinear global instability of a fluid flow through a collapsible channel

TL;DR

This paper investigates the linear and weakly nonlinear global instability of fluid flow in a collapsible channel, revealing bifurcation behaviors and the effects of membrane tension and Reynolds number through analytical models.

Contribution

It introduces a variant of an existing one-dimensional model using a weighted residual approach and analyzes bifurcation phenomena with weakly nonlinear theory.

Findings

Steady bifurcating flows occur at low membrane tension and Reynolds numbers.

Flow becomes time-periodic at high tension and Reynolds numbers.

Bifurcation types include transcritical saddle node and either subcritical or supercritical Hopf.

Abstract

Interactions between an internal flow and wall deformation occur in many biological systems. Such interactions can involve a complex and rich dynamical behavior and a number of peculiarities which depend on the flow parameter range. The aim of this paper is to present a variant (obtained via a weighted residual approach) of the averaged one-dimensional model derived by Stewart et al. "Local and global instabilities of flow in a flexible-walled channel," Eur. J. Mech. B/Fluids 28, 541-557 (2009)]. The asymptotic expansions for small Reynolds numbers of these two models, compared to the exact solution obtained from the lubrication approach, reveal some quantitative difference, even at higher Reynolds numbers. Qualitatively, the two models give similar results at least at a linear level. It is shown that for relatively low membrane tension (T), there are distinct regions in the (T,R) parameter space where steady bifurcating flows may occur. These flows can also be observed at vanishingly small Reynolds numbers combined with relatively high membrane tension. At sufficiently high T and R, the bifurcating flow is rather time periodic. A weakly nonlinear analysis is then performed in both cases leading to the derivation of evolution equations for the amplitudes of the bifurcating flows. The amplitude equations show that the saddle node bifurcation has a transcritical character while the Hopf bifurcation is either subcritical or supercritical, depending both on the mode number and membrane tension.

Figures

Tables

Table 1. Table 1: Takens-Bogdanov points tabulated for successive modes for both the WIM and IM.

Takens-Bogdanov points $(T,R)$ for the WIM and IM
Mode Number ($n$)	$1$	$3$	$5$	$7$
WIM	$(0.2237,15.30)$	$(0.0187,10.23)$	$(0.0064,9.83)$	$(0.0032,9.73)$
IM	$(0.1766,18.68)$	$(0.0146,12.63)$	$(0.0050,12.18)$	$(0.0025,12.05)$

Equations

$$u_x+v_y=0,$$

$$R(u_t+u{u}_x+v{u}_y+p_x)=u_{yy},$$

$$p_y=0,$$

$$\displaystyle(u,v)=(0,0)\quad\text{at}\quad\left\{\begin{array}[]{l}y=0,\ -L_{1}<x<1+L_{2}\\ y=1,\ -L_{1}<x<0\ \text{and}\ 1<x<1+L_{2}\end{array},\right.$$

$$(u,v,p)=(0,h_t,p_e-T{h}_{xx})y=h(x,t),0<x<1,$$

$$p_e(x)=p_u-\frac{p_u}{L_0}(x+L_1).$$

$$p=p_u-\frac{12}{R}(x+L_1)-T{h}_{xx}$$

$$R(u_t+u{u}_x+v{u}_y-T{h}_{xxx})=12+u_{yy}$$

$$q_x+h_t=0,$$

$$u(x,y,t)=\frac{q(x,t)}{{\int }_0^{h(x,t)}{u}_0\left(\frac{y}{h}\right)\mathrm{d}y}{u}_0\left(\frac{y}{h}\right)+R{u}_1(x,y,t)+O(R^2),$$

$$R\left(q_{t}-\frac{6}{5}\frac{q^{2}}{h^{2}}h_{x}+\frac{12}{5}\frac{q}{h}q_{x}-Thh_{xxx}\right)=12\left(h-\frac{q}{h^{2}}\right)+Ru_{1y}\Big{|}_{0}^{h}.$$

$$R\int_{0}^{h}f\left(\frac{y}{h}\right)\left\{\tilde{u}_{0t}+\tilde{u}_{0}\tilde{u}_{0x}+\tilde{v}_{0}\tilde{u}_{0y}-Th_{xxx}\right\}\mathrm{d}y=12\int_{0}^{h}f\mathrm{d}y+\int_{0}^{h}uf_{yy}\mathrm{d}y+fu_{y}\Big{|}_{0}^{h},$$

$$R\left(\gamma q_t-\frac{6}{5}\alpha \frac{q^2}{h^2}{h}_x+\frac{12}{5}\beta \frac{q}{h}{q}_x-Th{h}_{xxx}\right)=12\left(h-\frac{q}{h^2}\right),$$

$$\gamma q_t=-p_x-\frac{12}{R}q.$$

$$\displaystyle p\Big{|}_{x=0}=p_{u}-\left({\gamma}q_{t}+\frac{12}{R}q\right)L_{1},$$

$$\displaystyle p\Big{|}_{x=1}=\left({\gamma}q_{t}+\frac{12}{R}q\right)L_{2}.$$

$$\begin{array}[]{lll}x=0:&h=1,&Th_{xx}=L_{1}\left(\displaystyle\frac{12}{R}(q-1)+{\gamma}q_{t}\right),\\ x=1:&h=1,&Th_{xx}=-L_{2}\left(\displaystyle\frac{12}{R}(q-1)+{\gamma}q_{t}\right).\end{array}$$

$$q_0=h^3,$$

$$q_1=R{h}^3\left\{\frac{T}{12}{h}_{xxx}+\left(\frac{3}{4}\gamma +\frac{1}{10}\alpha -\frac{3}{5}\beta \right)h^3{h}_x\right\}.$$

$$h_t+(q_0+R{q}_1{)}_x=0.$$

$$u_0=6y(h-y),$$

$$u_1=R\left\{\frac{3}{2}h{h}_x(h^{3}y-2h{y}^3+y^4)+\frac{T}{2}{h}_{xxx}(hy-y^2)\right\},$$

$$q_0=h^3,$$

$$q_1=R{h}^3\left(\frac{3}{10}{h}^3{h}_x+\frac{T}{12}{h}_{xxx}\right).$$

$$q_1=R{h}^3\left(\frac{1}{4}{h}^3{h}_x+\frac{T}{12}{h}_{xxx}\right).$$

$${\tilde{q}}_x+{\tilde{h}}_t=0,$$

$$R\left(\gamma {\tilde{q}}_t-\frac{6}{5}\alpha {\tilde{h}}_x+\frac{12}{5}\beta {\tilde{q}}_x-T{\tilde{h}}_{xxx}\right)=12(3\tilde{h}-\tilde{q}),$$

$$\displaystyle\begin{array}[]{lll}x=0:&\tilde{h}=0,&T\tilde{h}_{xx}=L_{1}\left(\displaystyle\frac{12}{R}\tilde{q}+{\gamma}\tilde{q}_{t}\right),\\ x=1:&\tilde{h}=0,&T\tilde{h}_{xx}=-L_{2}\left(\displaystyle\frac{12}{R}\tilde{q}+{\gamma}\tilde{q}_{t}\right).\end{array}$$

$${\mathcal{L}}_s\tilde{h}=\frac{12}{R}{q}_1^{(s)},$$

$$\displaystyle\begin{array}[]{lll}x=0:&\tilde{h}=0,&T\tilde{h}_{xx}=\displaystyle\frac{12}{R}L_{1}q_{1}^{(s)},\\ x=1:&\tilde{h}=0,&T\tilde{h}_{xx}=-\displaystyle\frac{12}{R}L_{2}q_{1}^{(s)},\end{array}$$

$${\mathcal{L}}_s=T{\partial }_x^3+\frac{6}{5}\alpha {\partial }_x+\frac{36}{R}.$$

$$\left|\begin{array}{cccc}\frac{1}{3}& 1& 1& 1\\ -\frac{12L_1}{RT}& {\lambda }^2& {\bar{\lambda }}^2& (\lambda +\bar{\lambda }{)}^2\\ \frac{1}{3}& e^{\lambda }& e^{\bar{\lambda }}& e^{-(\lambda +\bar{\lambda })}\\ \frac{12L_2}{RT}& {\lambda }^2{e}^{\lambda }& {\bar{\lambda }}^2{e}^{\bar{\lambda }}& (\lambda +\bar{\lambda }{)}^2{e}^{-(\lambda +\bar{\lambda })}\end{array}\right|=0.$$

$${\mathcal{L}}_{\omega }{Q}_1=0,$$

$$\displaystyle\begin{array}[]{lll}x=0:&Q_{1x}=0,&\left\{T\partial_{x}^{3}+i{\omega}L_{1}\left(\displaystyle\frac{12}{R_{c}}+i\gamma\omega\right)\right\}Q_{1}=0,\\ x=1:&Q_{1x}=0,&\left\{T\partial_{x}^{3}-i{\omega}L_{2}\left(\displaystyle\frac{12}{R_{c}}+i\gamma\omega\right)\right\}Q_{1}=0,\end{array}$$

$${\mathcal{L}}_{\omega }=T{\partial }_x^4+\frac{6}{5}\alpha {\partial }_x^2+\left(\frac{36}{R_c}+\frac{12}{5}i\beta \omega \right){\partial }_x+\left(\frac{12}{R_c}i\omega -\gamma {\omega }^2\right)$$

$$det[{\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3,{\mathbf{v}}_4]=0$$

Full text

**Linear and weakly nonlinear global instability of a fluid flow

through a collapsible channel**

• Amaouche M, Di Labbio G*†* Laboratory of Cardiovascular Fluid Dynamics, Concordia University, Montréal, QC, Canada, H3G 1M8

Interactions between an internal flow and wall deformation occur in many biological systems. Such interactions can involve a complex and rich dynamical behavior and a number of peculiarities which depend on the flow parameter range. The aim of this paper is to present a variant (obtained via a weighted residual approach) of the averaged one-dimensional model derived by Stewart et al. “Local and global instabilities of flow in a flexible-walled channel,” Eur. J. Mech. B/Fluids 28, 541-557 (2009)]. The asymptotic expansions for small Reynolds numbers of these two models, compared to the exact solution obtained from the lubrication approach, reveal some quantitative difference, even at higher Reynolds numbers. Qualitatively, the two models give similar results at least at a linear level. It is shown that for relatively low membrane tension ($T$), there are distinct regions in the $(T,R)$ parameter space where steady bifurcating flows may occur. These flows can also be observed at vanishingly small Reynolds numbers combined with relatively high membrane tension. At sufficiently high $T$ and $R$, the bifurcating flow is rather time periodic. A weakly nonlinear analysis is then performed in both cases leading to the derivation of evolution equations for the amplitudes of the bifurcating flows. The amplitude equations show that the saddle node bifurcation has a transcritical character while the Hopf bifurcation is either subcritical or supercritical, depending both on the mode number and membrane tension.

• Data pertaining to this article will be made available from the authors upon reasonable request. * The authors have no conflicts of interest to declare. * Please cite as: Amaouche, M., & Di Labbio, G. (2016). Linear and weakly nonlinear global instability of a fluid flow through a collapsible channel. Physics of Fluids, 28(4), 044106.

© 2016 American Institute of Physics Publishing. This manuscript version is made available under the AIP Publishing LLC Transfer of Copyright Agreement, more information regarding usage terms can be found at https://publishing.aip.org/resources/researchers/rights-and-permissions/permissions/. This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in “Amaouche, M., & Di Labbio, G. (2016). Linear and weakly nonlinear global instability of a fluid flow through a collapsible channel. Physics of Fluids, 28(4), 044106.” and may be found at https://doi.org/10.1063/1.4946758.

Introduction

We investigate linear and weakly nonlinear instability of a viscous and incompressible Newtonian fluid driven by a prescribed pressure gradient through a finite length plane channel for which a segment of one wall is a membrane under longitudinal tension. This hydroelastic system is the plane schematic representation of the Starling resistor, first introduced by \citeAPedley92, in which a segment of flexible tubing is mounted between two rigid tubes and enclosed in a pressure chamber. The latter allows the external pressure acting on the flexible part of the tube to be controlled independently of the internal fluid pressure. When the elastic tube is subjected to a sufficiently large compression, it buckles asymmetrically along part of its length, leading to an increase of its compliance, thereby becoming potentially capable to develop large amplitude self-excited oscillations when the internal flow rate exceeds a certain critical value Bertram (\APACyear2003).

The propensity of flexible tubes to exhibit significant deformation has motivated a large number of investigations because of its relevance in a wide variety of physiological phenomena. For reviews of pertinent physiological applications, we refer to \citeAGrotberg04. These include noise-generating instabilities such as Korotkoff sounds produced during blood pressure measurement in the compressed brachial artery, snoring during inhalation, and flow limitation in lung airways such as wheezing during expiration. Self-excited oscillations are responsible for the generation of speech in the human larynx and bird song in the avian syrinx. It is presently known that self-excited oscillations in the Starling resistor arise in multiple distinct modes whose interactions can lead to a complex and rich dynamical behavior of the subsequent flow. Substantial advances have been achieved in understanding the mechanisms underlying these global instabilities, which are intimately related to the boundary conditions at the upstream and downstream ends of the flexible membrane.

Studies about global instabilities in finite-length asymmetric channels typically model the compliant wall as a tensioned membrane and neglect properties such as wall inertia, wall elasticity, and bending stiffness. The pressure is typically held fixed far downstream and either the flow rate or the pressure may be prescribed upstream. In either case, the system may exhibit a wide variety of self-excited oscillations. \citeAJensen03 identified, using high Reynolds number asymptotics, a mechanism for the generation of long wavelength and high frequency mode $1$ sloshing oscillations when the membrane is under high tension. They provided an explanation of this instability mechanism in terms of the fluid’s energy balance and found that this high frequency mode requires the upstream pressure to be prescribed and the upstream rigid channel segment to be shorter than that downstream. In these conditions, the in and out motion of the membrane forces significant oscillatory volume flow in the rigid parts of the channel and so energy can be transferred from the inflow to the oscillations. This high frequency mode has also been captured by the integral model (IM) derived by \citeAStewart09 who found that the base flow exhibits static and oscillatory global instabilities and showed how self-excited oscillations can arise entirely from wave reflections from the rigid parts of the system. Hence, an unstable global mode is composed of interacting local modes, each of which can be stable, with instability being driven primarily by the boundary conditions. Locally unstable modes are naturally expected to also contribute to the occurrence of global instability. Analyzing three different systems, \citeADoare06 have proposed a method to quantify the contributions of local instability and of wave reflection in the global instability of one-dimensional systems.

While quantitative accuracy can be guaranteed only for sufficiently small Reynolds numbers by the IM, the predictions for larger Reynolds numbers were shown to capture many important features computed by direct numerical simulation of the basic equations and by sophistical asymptotic analysis Stewart \BOthers. (\APACyear2010). The model derived by \citeAStewart09 was also employed by \citeAXu13, Xu14 to analyze the low frequency mode $1$ oscillations in the case of a fixed inflow rate. Changing the coefficient of convective inertia, they showed that the model’s predictions are qualitatively insensitive to the velocity profile assumption. Their weakly nonlinear analyses show in particular more than one type of unsteady behavior that can occur in the neighborhood of a particular organizing point. The results indicate that two possible states of sustained oscillations can arise, depending on whether the length of the downstream rigid segment is of the same order or much greater than the membrane length. For a detailed review of one-dimensional models we refer to \citeAPedley15.

As emphasized by many authors, the IM is a tractable model which performs well from small to moderate Reynolds numbers. However, despite its success in the nonlinear regime, this model does not incorporate the effects of inertia through the parabolic velocity assumption. Corrections to the parabolic velocity profile in order to coherently account for the inertial effects are therefore required for an accurate prediction of the linear instability threshold, at least for small inertia. Our purpose here is to obtain such a coherent first order model which also performs well for moderate Reynolds numbers as the IM does.

This paper is organized as follows. In Section II, we follow \citeAStewart09 to present the physical problem and its governing equations. In Section III, we derive the weighted integral model (WIM) and compare its asymptotic expansion with that obtained from the IM and the exact lubrication solution. Linear and weakly nonlinear stability of the base flow are examined in Sections IV and V, respectively. Representative numerical results from the two models are then given and compared in Section VI. Finally, a conclusion is offered in Sec. VII.

Physical Problem and Basic Equations

The physical system consists of a Newtonian incompressible fluid with density $\rho^{*}$ and dynamic viscosity $\mu^{*}$ flowing in a slender planar channel of total length $L^{*}_{0}$ (the asterix is used to denote dimensional quantities). A segment of one wall of the channel is made of a thin elastic membrane, of length $L^{*}$ and longitudinal tension $T^{*}$, mounted between two rigid segments with lengths $L^{*}_{1}$ upstream and $L^{*}_{2}$ downstream (see Fig. 1 for a schematic of the problem). To simplify matters, the other properties of the membrane (inertia, damping, and bending stiffness) are neglected. When the membrane is flat, the channel has a uniform width $d^{*}_{0}$ which is assumed much smaller than $L^{*}$. This introduces a small perturbation parameter $\epsilon=d^{*}_{0}/L^{*}$ which allows for appreciable simplifications in the formulation. The flow is driven by a prescribed constant pressure gradient between the channel exit where the pressure is $p^{*}=0$ and the channel entrance where $p^{*}=p^{*}_{u}$. We introduce a Cartesian coordinate system with the $x^{*}$ axis lying on the rigid wall such that $x^{*}$ increases downstream and $x^{*}=0$ corresponds to the abscissa of the upstream end of the membrane. The $y^{*}$ axis is oriented inside the flow. The smallness of $\epsilon$ implies a separation of scales between $x^{*}$ and $y^{*}$ which in turn dictates a different treatment of the streamwise and cross-stream momentum equations as in boundary layer theory. Hence, making use of a suitable scaling which will be defined shortly and neglecting terms of $O(\epsilon^{2})$ and higher, the dimensionless boundary layer type equations read

[TABLE]

where $u$ and $v$ are the streamwise and cross-stream velocity components, $t$ is time, and $R=U^{*}_{0}{d^{*}_{0}}^{2}/L^{*}{\nu}^{*}$ is the Reynolds number; ${\nu}^{*}$ being the kinematic viscosity and $U^{*}_{0}=p^{*}_{u}{d^{*}_{0}}^{2}/12{\mu}^{*}L^{*}_{0}$ being the averaged velocity of the base flow. Equations (1)-(3) are rendered dimensionless by scaling $x^{*}$, $y^{*}$, $t^{*}$, $u^{*}$, $v^{*}$, and $p^{*}$ with $L^{*}$, $d^{*}_{0}$, $L^{*}/U^{*}_{0}$, $U^{*}_{0}$, $U^{*}_{0}d^{*}_{0}/L^{*}$, and ${\rho}^{*}{U^{*}_{0}}^{2}$, respectively. The flow is subject to the no-slip and no-penetration boundary conditions at the walls and to the normal stress balance which determines the shape $y=h(x,t)$ of the membrane, namely,

[TABLE]

with $T=d_{0}^{*}T^{*}/{\rho}^{*}{U^{*}_{0}}^{2}{L^{*}}^{2}$ being the dimensionless tension parameter and $p_{e}(x)$ a prescribed external pressure defined below. When the membrane is flat, one obtains the base flow with velocity $(u_{0},v_{0})=(6y(1-y),0)$ from the equilibrium of viscous forces and the prescribed constant pressure gradient. The internal pressure distribution is dependent only on $x$ owing to (3) and is given by $p(x)=p_{u}-p_{u}(x+L_{1})/L_{0}$ with $p_{u}/L_{0}=12/R$. For there to be no pressure difference on either side of the membrane when it is flat, the same pressure distribution is imposed outside, namely,

[TABLE]

Equation (2) subject to the dynamic boundary condition at the membrane combined with (8) provides the pressure distribution in the channel for $0<x<1$,

[TABLE]

The momentum equation then becomes (for $0<x<1$):

[TABLE]

Model

Integrating the continuity equation over the channel width and accounting for the kinematic condition on the membrane yields the global equation

[TABLE]

where $q$ is the local flow rate. This shows that the flow rate, like the shape of the membrane, is a relevant variable. Consequently, it appears quite natural to seek another equation with suitable boundary conditions to obtain a closed system for these two genuine variables. Such a system was first derived by \citeAStewart09 for the Starling resistor. The IM is obtained by performing direct integration of the momentum equation and using a parabolic velocity profile resulting from the balance between viscous forces and the prescribed pressure gradient. Assuming such a velocity profile amounts to neglecting the distortion induced by inertia. For small Reynolds numbers, a suitable velocity profile accounting for the inertia effects would be of the form

[TABLE]

where $u_{0}(y/h)=6(y/h)(1-y/h)$ and the correction $u_{1}$ is such that $\int_{0}^{h}u_{1}\mathrm{d}y=0$. In a calculation up to $O(R)$, the correction $u_{1}$ is not needed to express inertial terms which are already of $O(R)$, but must be taken into account in the viscous term for consistency. Hence, up to $O(R)$ the averaged momentum equation would read

[TABLE]

This is the IM derived by \citeAStewart09 complemented by the unknown quantity Ru_{1y}\Big{|}_{0}^{h}. The question then is what global procedure allows for the incorporation of $u_{1}$ without calculating it explicitly. This can be achieved by introducing a weighting function $f(y/h)$ before integrating the momentum equation. Hence, after two successive integrations by parts and using the no-slip condition at the walls, one obtains

[TABLE]

where the tilde refers to the leading order term of $u$ in (12). One can see that the viscous term can be rendered free from $u_{1}$ if the weighting function is chosen so that $f_{yy}=K(x,t)$ and $f(0)=f(h)=0$, $K(x,t)$ being an arbitrary function of $x$ and $t$. So $f$ is, apart from a multiplicative constant, nothing more than the base velocity $u_{0}(y/h)$.

Further, the function $K(x,t)$ is chosen in order to have the same viscous term as in Equation (13); i.e., $K(x,t)=-12/h^{2}$ and therefore $f=6(y/h)(1-y/h)$. The integrated weighted momentum equation then reads

[TABLE]

with $(\alpha,\beta,\gamma)=(1,1,1)$ for the IM and $(\alpha,\beta,\gamma)=(9/7,17/14,6/5)$ for the present model referred to as the WIM. It is to be noted (i) that the term $-12q/h^{2}$ arising from the weighted integration of the viscous force is exact in the sense that it is obtained without any assumption on the velocity profile, and (ii) up to $O(R)$, Equation (14) is the same as Equation (13) without the term Ru_{1y}\Big{|}_{0}^{h} which is related to the unknown correction $u_{1}$.

Now, following \citeAStewart09, we now derive the boundary conditions associated to Equation (14) when the pressures at the ends of the tube are prescribed. Where the channel has rigid walls, Equation (14) can be simplified since $h=1$, which yields $q_{x}=0$ and therefore $q=q(t)$. Then, one is led to

[TABLE]

Integrating this equation subject to the prescribed pressures at the entrance and the exit of the channel yields

[TABLE]

Combining Equations (16) and (17) with (9), one obtains the boundary conditions for the two variable system given by Equations (11) and (14),

[TABLE]

In order to give evidence for the difference in predicting criticality for small Reynolds numbers between the WIM and the IM, we compare their asymptotic expansion as $R\rightarrow 0$ with the asymptotic solution of the full problem given by Equations (1) and (2). Using Equation(14) with $q$ being expanded in asymptotic series in $R$ (i.e., $q=q_{0}+Rq_{1}+\ \cdots$) and equating like powers of $R$, one obtains for the first two orders

[TABLE]

A \citeABenney66 type equation for the membrane shape is then derived from the continuity equation as

[TABLE]

In the range of Reynolds numbers where the Benney equation applies, the dynamics of the flow is slaved to its kinematics; i.e., the flow rate is adiabatically slaved to the membrane shape and it varies with $x$ and $t$ only through the dependence of the membrane shape on these variables.

Performing a similar asymptotic expansion for the velocity in (1) and (2) yields

[TABLE]

which correspond to

[TABLE]

Setting $(\alpha,\beta,\gamma)=(9/7,17/14,6/5)$ in (20), we retrieve the same correction for $q_{1}$ as that given by (24). It is worth noting that the lubrication approach requires the explicit calculation of the first order correction of the velocity profile whereas the WIM does not. In spite of this, the two approaches yield the same results for the correction $q_{1}$, at least up to $O(R)$, while the IM provides a slightly different expression, namely,

[TABLE]

Linear Stability

Following common practice, we consider linear stability of the basic solution $(h,q)=(1,1)$ by setting $(h,q)=(1,1)+(\tilde{h},\tilde{q})$, where $(\tilde{h},\tilde{q})$ are small perturbations satisfying the following linear system:

[TABLE]

The linearity of the system given by Equations (25)-(29) permits solutions proportional to $e^{{\sigma}t+{\lambda}x}$. The fixed point $(\tilde{h},\tilde{q})=(0,0)$ may lose its stability either in favor of another fixed point via a saddle node bifurcation or in favor of a limit cycle via a Hopf bifurcation, depending on the parameter values. The former occurs when a real eigenvalue $\sigma$ crosses the imaginary axis ($\sigma=0$) and the second corresponds to a pair of complex conjugate eigenvalues crossing.

Steady-state bifurcation

In this case, a secondary fixed point bifurcates from the basic state. This corresponds to a steady flow governed by equations which can be derived from (25) to (29) by ignoring the time dependence. Consequently, $\tilde{q}$ is an arbitrary constant, say $q_{1}^{(s)}$, and $\tilde{h}$ satisfies

[TABLE]

where $\mathcal{L}_{s}$ is the linear steady operator given by

[TABLE]

The solution is then $\tilde{h}\equiv q_{1}^{(s)}h_{1}^{(s)}$, where $h_{1}^{(s)}=1/3+ae^{{\lambda}x}+\bar{a}e^{\bar{\lambda}x}+be^{-(\lambda+\bar{\lambda})x}$, $\lambda$ being one of the two complex roots of the characteristic equation and $a$ and $b$ are complex and real constants, respectively. Nontrivial solutions exist provided

[TABLE]

Eigenvalue relation (35) yields neutral stability curves in the $(T,R)$ plane.

Hopf bifurcation

Instead of a single eigenvalue being at the origin, we assume now that if $R=R_{c}(T)$, $\sigma=\pm i\omega$ ($\omega>0$), which means that a pair of complex conjugate eigenvalues crosses the imaginary axis, while others are still in the left half-plane, resulting in a bifurcation into $2\pi/\omega$ time-periodic solutions. Along the bifurcation curve, setting $\partial_{t}=i\omega$, we have $\tilde{h}\equiv h_{1}^{(o)}=AQ_{1x}(x)e^{i{\omega}t}+C.C.$ and $\tilde{q}\equiv q_{1}^{(o)}=-i{\omega}AQ_{1}(x)e^{i{\omega}t}+C.C.$ where $A$ is an arbitrary constant and $Q_{1}$ is an eigenfunction of the following eigenvalue problem:

[TABLE]

where $\mathcal{L}_{\omega}$ is the linear operator given by

[TABLE]

Hence, $Q_{1}(x)=\sum_{n=1}^{4}c_{n}e^{ik_{n}x}$ with the $k_{n}$ being the roots of the quartic characteristic equation corresponding to (36). In order for a nontrivial solution to exist, the determinant of the linear system formed from the boundary conditions (39) must vanish

[TABLE]

with

[TABLE]

This gives the curve $R=R_{c}(T)$ along which the fixed point $(0,0)$ loses its stability in favor of a limit cycle. Representative results on the stability boundaries in the $(T,R)$ parameter space are shown in Fig. 2, illustrating a noticeable discrepancy in the predictions of the stability boundaries between the WIM and the IM.

Weakly Nonlinear Analysis

Beyond criticality and far enough from particular points where more than one mode are simultaneously destabilized, the amplitude of the most unstable global mode grows to a size sufficient to excite nonlinear self-interactions which retard its exponential temporal growth until a new flow sets in.

Steady bifurcation flow

Here, the aim is to correct the linear description of the bifurcating steady flow by accounting for the self-interactions of the steady global mode which is known to oppose to its exponential growth. To do this, instead of assuming that the amplitude of the critical global mode remains arbitrarily small, we consider now that it evolves slowly with time. Hence, it appears natural to introduce a fast time $t$ and a slow time $\tau={\delta}t$ where $\delta\ll 1$ measures the smallness of the flow amplitude near criticality. Then expanding the perturbations $(\tilde{h},\tilde{q})$ of the base state as $\tilde{h}={\delta}q_{1}^{(s)}h_{1}^{(s)}+{\delta}^{2}h_{2}^{(s)}+\ \cdots$, $\tilde{q}={\delta}q_{1}^{(s)}+{\delta}^{2}q_{2}^{(s)}+\ \cdots$, and $R=R_{s}(1+{\delta}R_{1})$, substituting $\partial_{t}$ by $\partial_{t}+\delta\partial_{\tau}$ in Equation (26), setting $\partial_{t}\equiv 0$ and separating terms of like powers of $\delta$, one first obtains the linear problem (30) at $O(\delta)$ and the following problem for $(h_{2}^{(s)},q_{2}^{(s)})$ at $O({\delta}^{2})$:

[TABLE]

subject to

[TABLE]

Here, the amplitude of $q_{1}^{(s)}$ is assumed to depend on the slow time $\tau$. Thus, solving (41) for $q_{2}^{(s)}$ gives

[TABLE]

The integration constant is dropped since it can be incorporated in the definition of $q_{1}^{(s)}(\tau)$. Then substituting the above expression for $q_{2}^{(s)}$ in the right hand side of (42) yields

[TABLE]

where

[TABLE]

The kernel of the operator $\mathcal{L}_{s}$ is not empty, therefore, for a solution to (45) to exist, the compatibility condition $S(x)\in Image(\mathcal{L}_{s})$ is required. Practically, this condition is expressed by invoking the solvability theorem for boundary value problems; i.e., the orthogonality (in the sense of some inner product to be defined) of $S(x)$ and the spanning function $h^{(a)}$ of the kernel of the adjoint operator $\mathcal{L}^{(a)}_{s}$ associated to $\mathcal{L}_{s}$, namely,

[TABLE]

where the inner product is defined as (see Appendix A for derivation)

[TABLE]

with

[TABLE]

Equation (47) along with Equations (48), (49), and (43) yields

[TABLE]

Combining this equation with (46) yields the following amplitude equation:

[TABLE]

which reads, after the change ${\delta}q_{1}^{(s)}\rightarrow q_{1}^{(s)}$,

[TABLE]

This equation characterizes a transcritical bifurcation with stability exchange. The coefficient $\alpha_{s}(R/R_{s}-1)$ is the exponential growth rate of the critical monotonic global mode while $\beta_{s}$ measures its quadratic self-interaction which is known to moderate its exponential growth during the linear stage.

Oscillatory bifurcation flow

Here, we apply the procedure developed above to the Hopf bifurcation with frequency $\omega$ corresponding to $R=R_{c}(T)$. It will turn out in the following that slightly past the threshold, the relevant slow time is $\xi={\delta}^{2}t$ and the corresponding expansion of $R$ is $R=R_{c}(1+{\delta}^{2}R_{2}+\ \cdots)$. The slow time $\tau={\delta}t$ and the $O(\delta)$ term in the expansion of $R$ would cancel due to the absence, at $O({\delta}^{2})$, of resonant terms. Then substituting for $\partial_{t}$ the operator $\partial_{t}+{\delta}^{2}\partial_{\xi}$, expanding $(\tilde{h},\tilde{q})={\delta}(h_{1}^{(o)},q_{1}^{(o)})+{\delta}^{2}(h_{2}^{(o)},q_{2}^{(o)})+{\delta}^{3}(h_{3}^{(o)},q_{3}^{(o)})+\ \cdots$, and identifying successive powers of $\delta$ in Equation (14) yields a sequence of three linear problems. The former is the marginal problem for $(h_{1}^{(o)},q_{1}^{(o)})$ whose solution is already calculated. At $O(\delta^{2})$ we get

[TABLE]

with boundary conditions

[TABLE]

The nonhomogeneous term in (54) is related to the eigensolution ($h_{1}^{(o)},q_{1}^{(o)}$) obtained in Section IV.B, namely,

[TABLE]

We note that (54) does not contain time-resonant terms, therefore it can be solved without any subsidiary condition. The solution of (54) is decomposed into steady and oscillatory parts as

[TABLE]

Solving (53) and (54) for $H_{20}(x)$ and $Q_{20}$, one obtains

[TABLE]

with

[TABLE]

where

[TABLE]

and $Q_{20}$, $c_{21}$, $\bar{c}_{21}$, and $c_{22}$ are yet to be determined constants by forcing the solution to satisfy the boundary conditions. We are led to a fourth order linear nonhomogeneous algebraic system of the form

[TABLE]

with

[TABLE]

[TABLE]

Now ignoring the steady term in (53)-(55), one obtains for $Q_{22}$,

[TABLE]

subject to

[TABLE]

In Equation (63), $\mathcal{L}_{2\omega}$ denotes the operator obtained from $\mathcal{L}_{\omega}$ in Equation (40) once $\omega$ is replaced by $2\omega$. Expressing $Q_{22}(x)$ in the form

[TABLE]

where $cof_{4n}(\mathbf{W})$ refers to the cofactor of the matrix $\mathbf{W}$ below resulting from deletion of the $4$th row and $n$th column and $a_{2n}$ are yet to be determined constants by using the above boundary conditions,

[TABLE]

So we get the following system:

[TABLE]

where

[TABLE]

The fast time dependence of quadratic self-interactions of marginal oscillatory global modes never has the frequency $\pm\omega$ which would lead to resonance. Accordingly, no quadratic terms occur in the amplitude equations. Therefore, the evolution equation for the amplitude $A(\xi)$ has to be formed at the next order which reads

[TABLE]

where $N.R.T.$ stands for non-resonant terms. They do not enter in the following calculation, so we refrain from giving them here because they are too cumbersome. The resonant term has the form

[TABLE]

where $S_{31}$, $S_{32}$, and $S_{33}$ are expressed in Appendix B.

Solving the system (67) and (68) by considering only the resonant source term, we can eliminate $h_{3}^{(o)}$ from (67) by setting $\partial_{t}=i\omega$ and $q_{3}^{(o)}=-i{\omega}Q_{3}(x)e^{i{\omega}t}$, Equation (68) then reduces to

[TABLE]

subject to

[TABLE]

where

[TABLE]

[TABLE]

The evolution equation for $A(\xi)$ will be found once again by invoking the Fredholm solvability condition (see Appendix C). One obtains the so-called Landau equation

[TABLE]

which can be rewritten by the change ${\delta}A\rightarrow A$ as

[TABLE]

Here the coefficient $\beta_{c}$ measures the cubic self-interaction of the critical oscillatory global mode. Its sign indicates the subcritical or supercritical nature of the bifurcation.

Numerics

Following \citeAStewart09, we assume throughout this section that $L_{2}$ is much greater than $L_{1}$ and set $L_{2}=10$ and $L_{1}=1$. Besides the fact that this configuration promotes oscillatory instability, as shown by these authors, it also allows for comparisons to be made. The curves drawn in Fig. 2a compare the neutral stability boundaries (corresponding to a real eigenvalue crossing the imaginary axis) in the $(T,R)$ parameter space obtained from the WIM (continuous lines) and the IM (dashed lines). They exhibit a noticeable discrepancy in the prediction of the monotonic instability threshold beyond some level of inertia. The departure between the two models becomes asymptotically small as $R\rightarrow 0$. Shaded regions indicate where the base flow is unstable against monotonic (static) modes and where a bifurcating steady flowsets in via a transcritical bifurcation. For small to moderate values of $T$, the onset of instability appears by decreasing or increasing $R$ depending on the shape (mode number) of the unstable mode. Each unstable region is bounded by two neutral stability branches termed “$n$stat” and “$(n+1)$stat” ($n=1,3,5,7,...$) indicating the mode number of the corresponding unstable mode. At the intersection points (marked by crosses), the modes $n$ and $n+1$ coalesce. The mode number $n$, which naturally decreases with increasing $T$, refers to the number of extrema in the spatial structures of the critical modes as shown in Fig. 5a. It is to be noted that monotonic instability occurs primarily for relatively small values of $T$, except when $R\rightarrow 0$ where higher values of $T$ are required to promote instability. A bifurcating steady flow may then take place in narrow bands provided the Reynolds number exceeds critical values which are approximately $2$ and $3$ for the modes $3$, …, $8$. Outside these bands, the pressure exerted by the tensioned membrane overcomes the destabilizing effect of inertia. However, for high enough Reynolds numbers inertia dominates the pressure effects and the instability changes character to become oscillatory and remains so for all values of $T$ and $R$.

Besides the stationary bifurcating flow, a time-periodic flow may also bifurcate from the base solution via a Hopf bifurcation. In the $(T,R)$ plane, this occurs along marginal curves termed “$n$osc” (see Fig. 2b), $n$ being the oscillatory mode number which also refers to the number of extrema in the shape of the critical mode (see Fig. 5b). Fig. 2b shows a strong dependence of the marginal Reynolds number $R_{c}$ on the membrane tension as $T$ falls to small values whereas it remains almost constant as $T$ is increased beyond some limit of order unity, regardless of the mode number. The base flow loses stability, in favor of a $2\pi/\omega$ time-periodic flow, when the critical branches are crossed by increasing $R$ with $T$ being fixed. Crosses indicate points where two critical modes of mode numbers ($2$,$3$), ($4$,$5$), and ($6$,$7$) coalesce. For example, modes ($2$,$3$) coalesce at $(T,R)=(0.03198,4.8388)$ or $(0.02493,6.1429)$ when using the WIM or IM, respectively. Where oscillatory and static critical curves intersect each other, monotonic and oscillatory modes interact and both bifurcate, which may result in complex dynamics and transitions. These co-dimension $2$ points are denoted by open and dark circles in Fig. 2. Points represented by dark circles mark the beginning of the marginal branches where the Hopf frequency vanishes (double-zero eigenvalue). Their proximity complicates the overall bifurcation structure which can be described by the Takens-Bogdanov normal form Guckenheimer \BBA Holmes (\APACyear1983). Table 1 compares, for the sake of illustration, the coordinates of the Takens-Bogdanov points in the $(T,R)$ plane obtained from the WIM and IM. Highly elaborate local and weakly nonlinear analyses by \citeAXu13 reveal how intricate the dynamics is around the particular Takens-Bogdanov point ($R\rightarrow\infty$, $T=T_{20}$), where oscillatory mode $2$ emerges through an interaction between two static eigenmodes when $L_{2}=O(1)$, before ultimately growing to large amplitude. The same organizing center was also investigated by \citeAXu14 but in the case of a longer downstream channel. This additional degree of freedom enables the system to sustain oscillations arising from resonance between mode $2$ and mode $1$ unstable oscillations.

Marginal curves for oscillatory modes $1$, $2$, $4$, and $6$ asymptote to $R=13.07$, $79.61$, $172.17$, and $308.76$, respectively, as $T\rightarrow\infty$. These values are found, from an asymptotic analysis similar to that employed in \citeAStewart09, to be $\gamma^{-1}$ times the corresponding values obtained by the IM. The variations of marginal mode frequencies $\omega$ with $T$ are plotted in Fig. 3. Rapid variations which are first observed for small $T$ are followed by asymptotic behavior of the form $a_{n}(\gamma^{-1}T)^{1/2}$, where $a_{n}$ indicates the sequence of eigenvalues of the leading order problem associated to (36) and (39) for sufficiently high values of $T$ (for the first two modes, one has $a_{1}=3.1807$ and $a_{2}=24.1251$ for example), namely,

[TABLE]

In Fig. 4, we compare the predictions of the critical Reynolds number and the corresponding critical frequency of the WIM and the IM with those given by direct numerical simulation in the same conditions explored by \citeAStewart10 (with the present scaling, the corresponding rigid channel lengths are $L_{1}=0.5$ and $L_{2}=3$). It is clear that the model of \citeAJensen03 is the most appropriate one to predict the instability threshold for the mode $1$ oscillation at sufficiently high membrane tension. At lower tension, less than about $100$, the IM and WIM are rather better candidates. The WIM is expected, by virtue of its derivation, to perform, if not better, at least as well as the IM does. However, even if good and equivalent estimates of the critical frequency are predicted by the two one-dimensional models, we note that the IM predicts the critical Reynolds number more accurately than the WIM. This apparent loss of accuracy of the WIM is expected to be caused by the fact that the viscous term $u_{xx}$ is incorporated in the DNS but is not taken into account in the models. Incorporating this second order term and using a suitable velocity profile for coherence would be a possible way to better appreciate the relative accuracies of the two integral models.

Typical monotonic and oscillatory eigenfunctions corresponding to modes $1$ through $4$ are represented in Fig. 5. As emphasized by \citeAHeil03, we observe that the steady membrane motions have larger amplitudes towards the downstream end of the membrane while time-periodic flows yield larger flow rates at the upstream end.

Fig. 6a depicts the coefficient $\alpha_{s}$ (the growth rate scaled by $R/R_{s}-1$) of the critical monotonic modes, as a function of $T$; it is important to keep in mind that instability corresponds to positive or negative values of $\alpha_{s}$, depending on the mode number. It is noticed that only modes with even mode number display a maximum whose value decreases with $T$. Close to Takens-Bogdanov points, the growth rate exhibits a singular behavior and changes sign on either side of those points. This a priori has no physical sense since $\alpha_{s}$ must remain negative everywhere on the branch $1$stat for instability to occur. The coefficient $\alpha_{s}$ must again be negative on almost the entirety of the branches corresponding to critical modes of odd number. This anomaly is expected to be related to the complexity of the dynamics arising from interaction of static and oscillatory modes near Takens-Bogdanov points. A better understanding of the physics around such points requires a local analysis similar to that developed by \citeAXu13 in the high $R$ regime. This study is in due course but is beyond the scope of the present paper. Plots in Fig. 6b represent the variations with $T$ of the eigenvalue $q_{1l}^{(s)}$ scaled by $(R/R_{s}-1)$, it is proportional to the amplitude of the static membrane deformation when the equilibrium state is reached. Values of $T$ are limited only, as mentioned above, to those leading to positive values of growth rate $\alpha_{s}(R/R_{s}-1)$. For those values of $T$, the product $(R/R_{s}-1)^{-1}q_{1l}^{(s)}$ is always negative which means that the static divergence increases or decreases the mean flow rate, depending on the marginal mode number; it decreases rapidly with $T$ except for the mode $2$stat for which the rate of decrease is slower for sufficiently high values of $T$.

Fig. 7a depicts the growth rate $\alpha_{c}$, scaled by ($R/R_{c}-1$), of the oscillatory modes against $T$. It can be observed that the degree of excitation follows the mode number; i.e., the higher the mode number, the larger the excitation, especially close to the Takens-Bogdanov points. Far from these points, the growth rate falls rapidly as $T$ increases for all modes except the mode $1$osc whose excitation remains almost constant (it is insensitive to variations of $T$). For this reason, the mode $1$osc is the dominant mode for values of $T$ greater than $0.2237$. The equilibrium amplitudes $A_{l}$ of the bifurcating oscillating modes, factored by $(R/R_{c}-1)^{-1}$, are plotted versus $T$ in Fig. 7b. Of utmost importance is the fact that the amplitudes of all the oscillating modes become singular at one or more particular values of $T$, namely, $T=2.9530$ for $1$osc and $T=0.0195$, $0.0291$, $0.0496$, $0.1052$, and $6.2218$ for $2$osc, for example. This breakdown is associated to the vanishing of the coefficients $\beta_{c}$, accounting for the cubic interaction of the modes when $R$ is slightly above $R_{c}$. This can be understood as a defect of the amplitude equations that might be alleviated by improving it through incorporating quintic order terms in view of the symmetry $A\rightarrow-A$. It is also worth noting that the product $(R/R_{c}-1)^{-1}|A_{l}|^{2}$ is negative in some ranges of $T$. For these values of $T$, the cubic term does not provide saturation of the linear exponential growth. As emphasized above, quintic order terms should be added to complete the amplitude equations, leading then to subcritical bifurcations with hysteresis phenomena which is associated with bistability. Stable time-periodic solutions also exist in some range of $R$ below $R_{c}$ and can be reached either via finite amplitude perturbations or by starting from a time-periodic regime and decreasing $R$ below $R_{c}$ for a fixed $T$. The subcritical character of the Hopf bifurcations beyond some limit of $T$ is in accordance with the asymptotic result obtained by \citeAStewart09.

From a quantitative point of view, performing a similar expansion to the one developed in Appendix A of \citeAStewart09, we found that the increase in the mean flowrate is $\beta$ times the one provided by the IM and the amplitude of the critical oscillatory mode is related to the one provided by the IM by a factor $\gamma^{-1/2}$. Finally, it is to be noted that at the linear and weakly nonlinear levels the differences between the WIM and the IM are only of quantitative character.

Summary

Linear and weakly nonlinear stability of planar Poiseuille flow in a channel with a flexible membrane as part of one of its walls is examined using a variant of the integral model. Despite the success of the IM in the nonlinear regime, it does not accurately predict neutral stability boundaries and introduces an error, typically of the order $20$% (due primarily to the presence of the coefficient $\gamma$ in the model equations) in the critical values of $T$ and $R$. The present model corrects this discrepancy but does not introduce any qualitative change in the IM predictions at the linear level. A weakly nonlinear analysis has revealed the nature of the two bifurcations encountered in the explored parameter range. Hopf bifurcations are subcritical for sufficiently large values of $T$ for all oscillating modes; this is in conformity with the asymptotic calculation as $T\rightarrow\infty$. For values of $T$ smaller than about $6$, Hopf bifurcations are supercritical or subcritical depending on both $T$ and the mode number.

Appendix A Adjoint Problem & Solvability Condition for (45) subject to (43)

The adjoint equation associated to Equation (45) can be derived by premultiplying it by some function $h^{(a)}$ and transforming the integral $\int_{0}^{1}h^{(a)}\cdot\mathcal{L}_{s}h_{2}^{(s)}\mathrm{d}x$ by repeated integration by parts and using the boundary condition $h_{2}^{(s)}(0)=h_{2}^{(s)}(1)=0$, one obtains

[TABLE]

where $\mathcal{L}^{(a)}_{s}$ is obtained from $\mathcal{L}_{s}$ by performing the change $\partial_{x}\rightarrow-\partial_{x}$. Now we define the adjoint operator $\mathcal{L}_{s}^{(a)}$ associated to the original problem (41)-(43) by $\mathcal{L}^{(a)}_{s}h^{(a)}=0$, $h_{x}^{(a)}(1)=h_{x}^{(a)}(0)=0$ which gives the adjoint eigenfunction

[TABLE]

where the boundary conditions $h_{x}^{(a)}(1)=h_{x}^{(a)}(0)=0$ imply

[TABLE]

and $c_{0}$ is a nonzero real constant which can be set equal to unity. Equation (A.1) then takes the reduced form

[TABLE]

which expresses the orthogonality of $h^{(a)}$ and the source terms in (42) and (43). Note that only two boundary conditions are sufficient to define the adjoint problem because $R_{s}$ is also an eigenvalue of $\mathcal{L}_{s}^{(a)}$ corresponding to $h^{(a)}$. It is possible to define the adjoint problem in another manner, namely by imposing a third boundary condition (a linear relation between $h^{(a)}(1)$ and $h^{(a)}(0)$, for example). This would lead to another spectrum of the adjoint operator and a different adjoint solution. A third boundary condition is however not required in the present case since Th_{2xx}\Big{|}_{0,1} is known at the boundaries.

Appendix B Expression of Source Terms in (69)

[TABLE]

Appendix C Adjoint Problem & Solvability Condition for (70) subject to (71)

In order to derive a solvability condition for the problem (70)-(71), we must define a suitable scalar product and the associated adjoint problem. We first define the adjoint operator of $\mathcal{L}_{\omega}$ as $\mathcal{L}_{\omega}^{(a)}=\bar{\mathcal{L}}_{\omega}^{*}$ with $\mathcal{L}_{\omega}^{*}$ being obtained from $\mathcal{L}_{\omega}$ by the change $\partial_{x}\rightarrow-\partial_{x}$. Let $Q_{a}$ be the eigenfunction of $\mathcal{L}_{\omega}^{(a)}$, i.e.,

[TABLE]

satisfying the following boundary conditions:

[TABLE]

where

[TABLE]

are the factors of the first $x$ derivative in $\mathcal{L}_{\omega}$. Then upon multiplying (70) by $\bar{Q}_{a}$ and performing successive integrations by parts yields the solvability condition

[TABLE]

which expresses the orthogonality of the source term $(S_{3}(x),B_{0},B_{1})$ and the complex conjugate of the spanning vector of the kernel of the adjoint operator $\mathcal{L}_{\omega}^{(a)}$.