Stability and the transition to turbulence in the flow through conduits with compliant walls

Abstract

The mechanisms of destabilisation of the flow through soft-walled channels/tubes are qualitatively different from those in rigid-walled conduits. The stability depends on two dimensionless parameters, the Reynolds number $(\rho _f V_f h_f/ \mu _f)$ and $\varSigma = (\rho _f G h_f^{2}/\mu _f^{2})$, where $\rho _f$ and $\mu _f$ are the fluid density and viscosity, $h_f$ and $V_f$ are the fluid length and velocity scale and $G$ is the wall elasticity modulus. There is an instability at zero Reynolds number when the dimensionless parameter $\varGamma = (\mu _f V_f / h_f G)$ exceeds a critical value. The low-Reynolds-number instability of the Couette flow past a compliant surface is well understood, and has been confirmed in experiments, but that in a pressure-driven flow is not completely understood. Two modes of instability at high Reynolds number have been predicted: the inviscid mode with an internal viscous layer, for which the transition Reynolds number scales as ${Re}_t \propto \varSigma ^{1/2}$; and the wall mode instability with a viscous layer at the wall, for which ${Re}_t \propto \varSigma ^{3/4}$. The wall mode instability has been observed in experiments at Reynolds number as low as 300 in a soft-walled tube and as low as 100 in a channel with one compliant wall, though the scaling of the transition Reynolds number differs from the theoretical prediction due to substantial wall deformation. Though the flow after transition shares many of the characteristics of hard-wall turbulence, it differs in significant ways, suggesting that soft-wall turbulence is a separate class distinct from hard-wall turbulence.

1. Introduction

Studies of the stability and transition in the flow past compliant surfaces have two distinct motivations. The first is drag reduction in marine and aerospace applications, which typically involve the external flow past large immersed bodies that are propelled through a fluid. In this case, the Reynolds numbers are very large, and the flow is turbulent over most of an object. The objective here is to either laminarise the flow or attenuate the turbulence by covering the object with a compliant ‘skin’, for efficient propulsion. The second motivation is transition in physiological flows in compliant conduits of dimension about 1 mm, where the Reynolds number is sufficiently low that the flow could be laminar. The central issue here is the effect of the wall compliance on the stability of the laminar flow, and on the nature of the flow after transition. In comparison to a laminar flow, the transition to turbulence could significantly enhance mixing in conduits of small dimension, where transport processes are often limited by slow molecular diffusion.

The present perspective considers the instability and the transition to turbulence in the flow through compliant tubes and channels. This field has benefited from an exchange of ideas with two adjacent areas of research: the external flows past objects with compliant surfaces and the flow in collapsible tubes and channels. Due to the difference in the length scales and the Reynolds numbers for external and internal flows, the methods of analysis have been different. Studies of external flows past compliant surfaces usually consider spatially developing boundary layer flows, and a spatial stability analysis is used to examine the growth of perturbations with downstream distance. In contrast, a fully developed unidirectional flow is usually considered in the flow through compliant conduits, and a temporal stability analysis is used to determine the transition Reynolds number. In the case of collapsible tubes, the focus has been on the shape oscillation of a compliant tube due to a difference between the internal and external pressures, and not on the laminar–turbulent transition itself. A brief summary of the stability and transition in external flows and collapsible tubes is first provided in this introduction, before proceeding to consider the stability of internal flows.

In marine and aerospace applications, the flow over large immersed objects is usually turbulent. Even when the flow is laminar at the upstream side, there is boundary layer separation, transition and turbulence downstream. It is desirable to reduce the drag force for more efficient propulsion, and there have been many studies carried out to examine whether it is possible to reduce drag by covering an object with a compliant coating. Experiments were conducted by Kramer (1960a,b, 1962) on dolphin-shaped objects, covered with viscoelastic materials with varying compliance, in tanks of water. These were based on the hypothesis that there could be additional dissipation of energy in the compliant material, resulting in the laminarisation of the flow around the object. The studies did find significant drag reduction of up to 40 % when a compliant surface was employed. The drag was found to first decrease and then increase as the dissipation in the wall material increased, suggesting that the properties of the surface could be tuned to attain maximum drag reduction. Later studies of other types of compliant surfaces (discussed in Carpenter & Garrad 1985) did not find a significant drag reduction, suggesting that this phenomenon is sensitive to the properties of the compliant surface. The promise of the use of compliant surfaces for drag reduction has not yet been realised in commercial applications. However, the initial observations led to a large number of theoretical studies in two distinct areas, the transition delay and turbulence attenuation due to compliant surfaces.

The linear stability of the flow past a compliant surface was first examined by Benjamin (1960, 1963) and Landahl (1962), who considered the flow past the ‘spring-backed wall’ model discussed in § 3.1. Here, the no-penetration condition applicable for a rigid surface is replaced by a relation between the pressure and the normal displacement of the wall. The no-slip condition is usually used for the tangential velocity at the wall, though there have been studies where the tangential displacement is related to the shear stress as well. Those authors identified three different modes of destabilisation in the flow past a compliant surface. The first is the class A mode which is a modification of the Tollmien–Schlichting mode for the flow past a rigid surface, modified due to wall compliance. These modes were found to be stabilised by wall compliance, but could be destabilised due to internal dissipation in the wall material. The class B modes are waves travelling at velocities close to those of free surface waves, and are considered a resonance effect. These are destabilised by an increase in wall compliance. The class C modes are analogous to the Kelvin–Helmholtz instability due to transfer of fluctuating energy from fluid to solid.

A subsequent study by Carpenter & Garrad (1985, 1986) classified these instabilities into two categories: the Tollmien–Schlichting mode which is present in the flow past a rigid surface and a class of instabilities collectively called flow-induced surface instability (FISI). The former is a modification of the instability in the flow past a rigid surface, and fluid viscosity is necessary for destabilising the flow. In contrast, the latter is a continuation of the inviscid modes at high Reynolds number, and fluid viscosity is not necessary for destabilising the flow. Carpenter & Garrad (1985, 1986) found that dissipation in the wall had a mild destabilising effect on the Tollmien–Schlichting mode, but stabilised the FISI. An important result was that a coalescence of the Tollmien–Schlichting mode and one of the FISIs could lead to a powerful new instability called the static divergence. Thus, these linear stability studies showed that although wall compliance could stabilise the Tollmien–Schlichting modes, it creates other instabilities that are not present in the flow past rigid surfaces.

A hydroelastic instability in the form of static divergence waves was observed in the experiments of Gad-el Hak, Blackwelder & Riley (1985), where a plate covered with a compliant surface was towed in a tank of water. Experiments were carried out for laminar, transitional and turbulent boundary layers. Static divergence waves were observed only for turbulent boundary layers when the velocity exceeded a transition value, and these were not observed for laminar boundary layers. This is in contrast to the experiments of Hansen & Hunston (1974, 1983) for a disk coated with a compliant material rotating in a tank of fluid, where static divergence waves were observed for both laminar and turbulent boundary layers. The static divergence waves on the compliant material caused flow modification similar to roughness elements on a rigid surface. The difference in the results of Gad-el Hak et al. (1985) and Hansen & Hunston (1974, 1983) has not been resolved so far, and more experimental work is required to make a connection between theoretical and experimental studies.

Studies of the flow through compliant tubes and channels can be broadly classified in two categories: the structural instability of a collapsible tube/channel and the flow instability in a conduit with compliant walls due to fluid–wall interaction. The present perspective is restricted to the latter. The structural instability in collapsible tubes is used to model phenomena such as venous collapse in the cardiovascular system and wheezing in the respiratory system. A typical configuration is the Starling resistor (Knowlton & Starling 1912), which is an elastic tube fixed between two rigid pipes with a constant pressure on the outside. The two-dimensional analogue of this was first studied by Pedley (1992). In its simplest manifestation, as the flow rate is increased for an inviscid flow, the pressure within the tube decreases and this could result in the collapse of the tube. However, a more detailed study reveals additional effects such as flow separation due to the downstream divergence in a collapsed tube (Cancelli & Pedley 1985). One- and two-dimensional models for the flow of varying complexity have been formulated (Grotberg & Jensen 2004), and these exhibit bifurcations between normal, buckled (Heil & Pedley 1996), distended or collapsed as well as self-excited oscillations (Jensen 1990) and high-frequency flutter (Gavriely et al. 1989). These collapsed states and oscillations have been carefully mapped out in experiments (Bertram, Raymond & Butcher 1989; Bertram, Raymond & Pedley 1990). The high-frequency flutter has been related to wheezing (Gavriely et al. 1989) and to the Korotkoff sounds that are used to make clinical diagnoses (Ur & Gordon 1970; Bertram, Raymond & Butcher 1989). The reader is referred to a review article (Grotberg & Jensen 2004) for further details on this subject.

In the case of collapsible tubes and channels, models are used to relate the internal pressure and the cross-sectionally averaged velocity, the latter being only a function of the streamwise coordinate. The mass conservation equation relates the time evolution of the area of cross-section to the cross-sectionally averaged velocity. The momentum conservation equation relates the substantial derivative of the cross-sectionally averaged velocity to the local pressure gradient, the viscous friction and other effects. The model relating the pressure difference across the wall of the tube and cross-sectional area includes the effect of tube elasticity and bending stiffness. Transitions, oscillations and bifurcations in the tube shape and flow velocity are then predicted using evolution equations for the cross-sectional area and the average velocity. The studies of collapsible tubes focus on the steady and periodic oscillations in the shape of the tube due to flows that may already be turbulent, and not on the transition from a laminar to a turbulent flow. In contrast to the study of stability in conduits with compliant walls, cross-sectionally averaged velocity and pressure fields are used in the flow through collapsible tubes, and the focus is on structural transitions and not on the laminar–turbulent transition.

It is instructive to discuss the instability in parallel flows through rigid channels and tubes to provide a basis for understanding the flow through compliant conduits. In experiments, the Reynolds number for the transition to turbulence is about 1200 (Patel & Head 1969) for a parabolic flow in a two-dimensional channel, and about 2100 (Reynolds 1883) for a cylindrical pipe. Since the transition takes place at high Reynolds number, it might naively be assumed that viscous effects can be neglected, and it is sufficient to consider the inviscid equations in order to predict transition. However, there are theorems such as the Rayleigh inflection point theorem which state that an inviscid flow can be unstable only if there is an inflection point somewhere in the flow. Since the parabolic laminar profile in a channel does not have an inflection point, an inviscid analysis predicts that the flow is always stable.

It turns out that the Tollmien–Schlichting instability for channel flow is due to the presence of an internal critical layer within the flow where viscous effects are important. This internal critical layer is of thickness ${Re}^{-1/3}$ smaller than the channel width at the location where the flow velocity is equal to the velocity of the waves. Here, ${Re}$ is the Reynolds number. The Tollmien–Schlichting mode becomes linearly unstable at a Reynolds number of about 5772 in a two-dimensional channel (Orszag 1971; Drazin & Reid 1981). In experiments, the transition to turbulence is observed at a Reynolds number of about 1200. This discrepancy is considered to be due to the highly subcritical nature of the instability – even though flow is unstable to infinitesimal perturbations at a Reynolds number of 5772, small but finite-amplitude perturbations become unstable at a much lower Reynolds number.

Stability analyses predict that the laminar flow in a cylindrical pipe is always stable. In experiments, the transition is observed at a Reynolds number of about 2100. This discrepancy is attributed to the highly subcritical nature of the bifurcation – even though the flow is stable to infinitesimal disturbances, it is unstable to finite-amplitude perturbations when the Reynolds number exceeds 2100 (Darbyshire & Mullin 1995; Draad, Kuiken & Nieuwstadt 1998; Hof, Juel & Mullin 2003; Mullin 2011). This hypothesis is supported by experiments which show that the flow could be maintained in the laminar state at Reynolds numbers up to $10^{5}$ if care is taken to avoid disturbances (Pfenniger 1961).

The nature of transition in channel and pipe flows is still poorly understood, over a hundred years after it was discovered. In experiments, the onset of transition at Reynolds number of about 2100 occurs due to localised disturbances, turbulent puffs and slugs, which originate at the wall of the channel of the pipe (Wygnanski & Champagne 1973). Transition due to disturbances imposed by injection and suction at the wall occurs at a higher Reynolds number in ‘clean flows’ without disturbances (Draad et al. 1998). Two broad theoretical explanations have been explored in the recent past. The first is the rapid algebraic growth of perturbations which are linearly stable due to the non-normal nature of the differential operator in the linear stability problem (Gustavsson 1991; Henningson, Lundbladh & Johansson 1993; Grossmann 2000). Although the algebraic growth of an isolated perturbation is followed by exponential decay much later on, it is assumed that there is onset of nonlinear interactions prior to decay, thereby sustaining turbulence.

Ideas from dynamical systems theory have been used to describe turbulence since the early attempts of Landau (1944) and Hopf (1948), who proposed that the flow undergoes a series of bifurcations after transition leading to a turbulent flow. This view has not been validated by later studies, and it is now accepted that the transition is discontinuous. Methods have been developed to identify low-dimensional models for turbulent flow (Berkooz, Holmes & Lumley 1993) using proper orthogonal decomposition of experimental/computed turbulent flows. For wall-bounded turbulent flows, the mechanism of turbulence generation is the roll-up of spanwise vortices into hairpin eddies and subsequent bursting of these eddies. Nonlinear and non-turbulent time-periodic travelling wave solutions of the Navier–Stokes equations have been identified in simulations (Nagata 1990; Waleffe 2001; Kerswell 2005; Kawahara, Uhlmann & van Veen 2012). Attempts have been made to relate the transitions among these states to the turbulence generation mechanism (Eckhardt et al. 2007). These solutions, labelled ‘exact coherent states’ (Graham & Floryan 2021), are not attractors, but are saddle points in a high-dimensional state space. The proposal is that the system travels close to these saddle nodes for relatively long periods of time, and therefore accurate identification of these solutions could capture the essential features of turbulent flows. However, a large number of these states have been identified for turbulent flows (Graham & Floryan 2021), and the total number of such states is not known. Moreover, the relation between these states and the transition to turbulence is not clear, since no state has yet been uniquely identified as the pathway to turbulence.

Linear stability analyses do predict instabilities in channels and tubes with compliant walls. Some of these are not just continuations of the instability in a rigid channel/tube, but are qualitatively different. They involve a dynamical interaction between the flow and the wall dynamics. In these cases, the transport of energy from the mean flow to the fluctuations takes place due to the shear work done at the interface, and not due to the Reynolds stresses within the fluid. In contrast to the situation for the flow in rigid conduits, experimental results for the critical Reynolds number for transition to a non-laminar state are consistent with the predictions of linear stability analysis. Thus, available evidence suggests that the transition observed in compliant channels/tubes is due to a linear instability to infinitesimal perturbations. In this sense, the transition in compliant conduits seems to be better understood in comparison with that in rigid conduits.

There have been very few experiments on transition and turbulence in internal flows preceding the theoretical studies that showed the existence of a linear instability. The reason could be that transition/turbulence was not considered important in biological flows at small dimensions and low velocities. The pioneering experiments in this area were carried out by Lahav, Eliezer & Silberberg (1973) and Krindel & Silberberg (1979) who first raised the possibility that, in comparison with rigid conduits, the transition could occur at a lower Reynolds number in conduits with compliant walls. It was only in the 1990s that linear stability studies identified a mechanism of destabilisation, the transport of energy from the mean flow to the fluctuations due to the shear work done at the interface, which is qualitatively different from that in rigid conduits. This led to further experiments with different geometries that are described in § 8.

The linear stability equations for a Newtonian fluid are briefly introduced in § 2 in order to clarify the notation and the scalings used. Along with a discussion of different wall models, a reasonably detailed description of the equations for a viscoelastic solid is provided in § 3. This is because the stability characteristics are known to be very sensitive to the solid equations at low Reynolds number, and erroneous stability results are obtained if care is not taken with regard to the many intricacies in the solid equations and boundary conditions. The Lagrangian and Eulerian formulations for a solid continuum are presented, and then the constitutive relations for a neo-Hookean solid are derived for the Lagrangian formulation. The stability analysis in the viscous limit at low Reynolds number is the subject of § 5. The different destabilisation mechanisms at high Reynolds number are discussed in § 6. Here, the asymptotic results for the inviscid and wall modes are first explained, and then the numerical results are discussed. Many weakly nonlinear stability analyses have been undertaken to understand the nature of the bifurcation after transition. The weakly nonlinear analysis is algebraically complicated, and so it is not discussed here, but the reader is referred to Drazin & Reid (1981) for an explanation of this procedure. The results of the weakly nonlinear analyses are discussed along with those of linear stability studies where appropriate. The experimental results are reviewed, and compared with the theoretical results, in § 8. A brief summary of the relatively small number of studies of simulations of turbulent flows is provided in § 9.

2. Fluid equations

Most studies are carried out with Newtonian fluids, where the governing equations are the incompressible Navier–Stokes mass and momentum equations:

(2.1)$$\begin{gather} \boldsymbol{\nabla} {\boldsymbol \cdot} {\boldsymbol v} = 0, \end{gather}$$

(2.2)$$\begin{gather}\rho_f \left( \frac{\partial {\boldsymbol v}}{\partial t} + {\boldsymbol v} {\boldsymbol \cdot} \boldsymbol{\nabla} {\boldsymbol v} \right) ={-} \boldsymbol{\nabla} p_f + \mu_f \nabla^{2} {\boldsymbol v}, \end{gather}$$

where ${\boldsymbol v}$ is the velocity, $p_f$ is the pressure, $\rho$ is the density and $\mu _f$ is the fluid viscosity. The stress tensor in the fluid is

(2.3)\begin{equation} \boldsymbol{\sigma}^{f} = \mu_f (\boldsymbol{\nabla} {\boldsymbol v} + (\boldsymbol{\nabla} {\boldsymbol v})^{T}).\end{equation}

The base state is usually a unidirectional flow, where the streamwise velocity $v_x$ is only a function of the cross-stream $(y)$ direction. Without loss of generality, the compliant wall is placed at the location $y=0$, and the other wall at the location $y = h$ could be a moving wall (for a Couette flow) or a stationary wall (for a Poiseuille flow). The mean velocity profile $\bar {v}_x$ for a fully developed unidirectional flow in a channel is

(2.4)\begin{equation} \bar{v}_x ={-} \frac{h^{2}}{4 \mu_f} \frac{\textrm{d} \bar{p}_f}{\textrm{d} x} \left( \frac{y}{h} - \frac{y^{2}}{h^{2}} \right) + \frac{V_w y}{h}, \end{equation}

where $(\textrm {d} \bar {p}_f / {\textrm {d} x})$ is the mean pressure gradient and $V_w$ is the velocity of the wall at $y=h$. The strain rate at the wall is

(2.5)\begin{equation} \dot{\gamma}_w ={-} \frac{h}{4 \mu_f} \frac{\textrm{d} \bar{p}_f}{\textrm{d} x} + \frac{V_w}{h}. \end{equation}

Though (2.4) is the exact solution for the steady and fully developed velocity profile, other approximate solutions of the mean velocity profile have also been used, along with the parallel flow approximation for flows that are slowly developing in the streamwise direction. This approximation has been used in Shankar & Kumaran (1999) for the stability of the developing flow in a cylindrical tube surrounded by a compliant wall and the stability of a converging flow in a tube with small angle of inclination. In the case of a developing flow, it is assumed that the wavelength of the perturbations is smaller than the length for flow development, so that the flow can be considered to be locally invariant in the flow direction. The ratio of the flow development length and the channel height or tube diameter is ${\sim }0.03 {Re}$, where $Re$ is the Reynolds number. Therefore, the steady parallel flow approximation is valid only for high Reynolds number, in the limit where the wavelength of perturbations is much smaller than the development length.

The flow in a converging/diverging channel is approximated as a parallel flow when the angle of inclination of the wall is small. For a rectangular channel or a cylindrical tube, the mean velocity profile (2.4) is determined by balancing the pressure gradient and the divergence of the viscous stress in the steamwise momentum equation, because the inertial term is identically zero for a steady fully developed unidirectional flow. For a slowly converging or diverging channel/tube with angle of wall inclination $\alpha \ll 1$, the correction to the viscous term is $O(\alpha )$, whereas the correction to the inertial terms is $O({Re} \alpha )$. Therefore, there could be a substantial modification of the mean velocity profile even when the angle of inclination is small, provided ${Re} \alpha \sim 1$.

In the linear stability analysis, perturbations are imposed on the base state in the form of plane waves in the streamwise ($x$) and spanwise ($z$) directions for the flow in a channel:

(2.6)\begin{align} \left. \begin{array}{c@{}} v_i = \bar{v}_i(y) + \tilde{v}_i(y) \exp{(\imath k_x x + \imath k_z z + s t)} = \bar{v}_i(y) + \tilde{v}_i \exp{(\imath k_x (x - c t) + k_z z)}, \\ p_f = \bar{p}_f + \tilde{p}(y) \exp{(\imath k_x x + \imath k_z z + s t)} = \bar{p}_f(y) + \tilde{p}_f \exp{(\imath k_x (x - c t) + k_z z)}, \end{array} \right\} \end{align}

where $\imath = \sqrt {-1}$, $k_x$ and $k_z$ are the wavenumbers in the streamwise and spanwise directions and $s$ is the growth rate of the perturbations. The flow is stable if the real part of $s$ is negative, indicating that perturbations decay exponentially, and unstable if the real part of $s$ is positive, indicating that perturbations grow exponentially. The growth rate $s$ is also written as $s = - \imath k_x c$, where $c$ is the wave speed. In this case, perturbations are stable if the imaginary part of $c$ is negative and unstable if the imaginary part of $c$ is positive.

These perturbations are substituted into the mass and momentum equations, which are then linearised in the perturbation amplitudes to obtain the linear stability equations:

(2.7)$$\begin{gather} \imath k_x \tilde{v}_x + \frac{\partial \tilde{v}_y}{\partial y} + \imath k_z \tilde{v}_z = 0, \end{gather}$$

(2.8)$$\begin{gather}\rho_f (s + \imath k_x \bar{v}(y)) \tilde{v}_x + \tilde{v}_y \frac{\textrm{d} \bar{v}}{\textrm{d} y} ={-} \imath k_x \tilde{p}_f + \mu_f \left(\frac{\textrm{d}^{2}}{{\textrm{d} y}^{2}} - k_x^{2} - k_z^{2} \right) \tilde{v}_x, \end{gather}$$

(2.9)$$\begin{gather}\rho_f (s + \imath k_x \bar{v}(y)) \tilde{v}_y ={-} \frac{\textrm{d} \tilde{p}_f}{\textrm{d} y} + \mu_f \left(\frac{\textrm{d}^{2}}{{\textrm{d} y}^{2}} - k_x^{2} - k_z^{2} \right) \tilde{v}_y, \end{gather}$$

(2.10)$$\begin{gather}\rho_f (s + \imath k_x \bar{v}(y)) \tilde{v}_z ={-} \imath k_z \tilde{p} + \mu_f \left(\frac{\textrm{d}^{2}}{{\textrm{d} y}^{2}} - k_x^{2} - k_z^{2} \right) \tilde{v}_z. \end{gather}$$

There have been relatively fewer studies of the stability of non-Newtonian fluids in conduits with compliant walls (Shankar & Kumar 2004; Chokshi & Kumaran 2007; Chokshi, Bhade & Kumaran 2015; Giribabu & Shankar 2017; Patne & Shankar 2019b). The models used and the salient results are described in § 7.

3. Wall models

The conduits are considered to have uniform cross-section in the base state, where the solid displacement field is considered to be unidirectional, steady and fully developed in most cases. It is shown in § 3.3.3 that for a steady fully developed flow, the displacement field in a viscoelastic solid wall is also fully developed and invariant in the streamwise direction (Gaurav & Shankar 2010). A similar result can be obtained for the flow in a tube (Gaurav & Shankar 2009). Therefore, it is valid to assume a base state with properties invariant in the flow direction. However, for other types of walls such as the spring-backed wall discussed in § 3.1, the wall displacement in the base state is not invariant in the flow direction. A stability analysis that assumes a channel/tube of infinite extent is an approximation that can be used only at high Reynolds number, where the inertial contributions to the stress perturbations are dominant.

3.1. Spring-backed plate

The spring-backed plate is a model that has been used extensively in the study of external flows past compliant surfaces (Benjamin 1960, 1963; Landahl 1962; Carpenter & Garrad 1985, 1986; Yeo & Dowling 1987), and less often for internal flows (Gajjar & Sibanda 1996; Kumaran 1996; Davies & Carpenter 1997a; Larose & Grotberg 1997; Shankar & Kumaran 2000; Thaokar, Shankar & Kumaran 2001). As the name suggests, the model consists of a deformable plate in contact with a fluid above, which is supported below by an array of springs and a fluid which acts as a damping element, as shown in figure 1(a). The surface is flat in the base state, and the deformation of the surface due to pressure fluctuations is described by the displacement field $u$ which is a function of the position along the surface. The equation for the displacement field is of the form

(3.1)\begin{equation} I \frac{\textrm{d}^{2} u}{\textrm{d} t^{2}} ={-} E u - D \frac{\partial u}{\partial t} + T \nabla_s^{2} u \pm p_f, \end{equation}

where $p_f$ is the difference in the fluid pressure on the surface and the pressure in the fluid supporting the surface from below; the positive sign on the right-hand side of (3.1) is applicable at the upper surface $y=y_u$ and the negative sign is applicable at the lower surface at $y=y_l$ in figure 1(b). Equation (3.1) is Newton's second law applied to a unit area of the surface, where $E$ is the spring constant, $I$ and $D$ are the effective inertia and damping coefficient per unit area, $T$ is the surface tension and $\nabla _s$ is the gradient operator along the undeformed surface. The no-slip boundary condition is used for the tangential velocity at the surface in the study of external flows past compliant surfaces.

Figure 1. (a) Schematic of a Couette flow between a rigid surface and a spring-backed wall and (b) pressure-driven flow between two spring-backed walls showing the deformable plate, base plate, spring and damping elements, the coordinate system and the displacement $u$.

The spring-backed plate model provides a simple relationship between the fluid pressure and the displacement of the surface. It can only be applied to high-Reynolds-number flows modelled using the inviscid approximation, since there is no tangential stress balance condition along the surface. For pressure-driven flows, there is a variation in the fluid pressure along the streamwise direction, which results in a variation in the displacement at steady state. Since this variation is neglected when the configuration is approximated as a plane channel or a cylindrical pipe, the spring-backed plate model can be used only when the surface displacement is much smaller than the cross-stream distance. The tangential stress at the surface also results in a tension in the interface along the tangential direction, which has to be balanced by a tangential spring force. Instead of a zero tangential velocity condition, some studies have used a wall model with a tangential displacement field, $u_T$, which is described by an equation of the form

(3.2)\begin{equation} I_T \frac{\partial^{2} u_T}{\partial t^{2}} ={-} E_T u_T - D_T \frac{\partial u_T}{\partial t} + \sigma_T, \end{equation}

where $\sigma _T$ is the shear stress, the subscript $T$ refers to the tangential direction and $E_T$, $D_T$ and $I_T$ are the spring, damping and inertia coefficients in the tangential direction. This type of interface condition was postulated in Thaokar et al. (2001) for a viscous flow and in Larose & Grotberg (1997) for the high-Reynolds-number flow in a compliant channel. Thaokar et al. (2001) found that inclusion of the tangential displacement in the model qualitatively alters the nature of the instability at zero Reynolds number.

The spring-backed plate model is useful because general requirements for the presence of unstable modes can be derived, such as the equivalents of Squire, Rayleigh and Fjørtoft theorems (Drazin & Reid 1981). An example of the derivation for a plane channel, along the lines of Yeo & Dowling (1987), is given in § 6.1. The equivalents for a pipe flow have been derived in Kumaran (1996) and Shankar & Kumaran (2000).

One important disadvantage of the spring-backed plate model with normal displacement, (3.1), is that the mean pressure gradient results in wall displacement that varies along the streamwise direction. The channel width or pipe diameter varies in the streamwise direction. This is typically not incorporated in the study of high-Reynolds-number flows, where it is assumed that the slope of the wall is small and the parallel flow assumption is valid. The advantage of using a model with a tangential displacement field, (3.2), is that the tangential spring force balances the wall shear stress locally, and the base state is a fully developed flow with a constant tangential displacement and zero normal displacement in the solid.

3.2. Membrane

Another model that has been used is the membrane model, where a tensioned membrane of infinitesimal thickness forms the wall of the channel, as shown in figure 2(a). A typical configuration considered in Kumaran & Srivatsan (1998) and Thaokar & Kumaran (2002) consists of a tensioned membrane stretched between two fluids, one or both of which are subject to a linear shear flow. In Stewart, Waters & Jensen (2009, 2010), there is only one fluid layer, and the pressure outside of the membrane is considered a constant. In this case, due to the shear stress exerted by the fluid on the surface, there is a variation in the membrane tension along the surface. In order to assume that the membrane properties are invariant along the streamwise direction in a normal mode analysis, it is necessary to consider a strongly pretensioned membrane in which the variation in the tension over a distance comparable to one wavelength is much smaller than the absolute tension. Kumaran & Srivatsan (1998) and Stewart et al. (2009, 2010) considered a relation similar to (3.1) for the membrane, without the inertia and spring forces, as the model for the membrane. The tangential velocity was set equal to zero at the membrane surface in this simple model.

Figure 2. The configuration and coordinate system for analysing the flow past a membrane (a) and the in-plane displacement field due to stress perturbations at the membrane surface (b).

A more sophisticated model was considered in Thaokar & Kumaran (2002), where the two-dimensional surface displacement field was defined along the surface of the membrane, as shown in figure 2(b). The base state is the tensioned state of the membrane with steady fluid flow. A material point at the location $(x_0, z_0)$ in the base state moves to a new location $(x,z)$ due to applied perturbations, and the two-dimensional displacement vector on the surface from the initial to the final location is designated ${\boldsymbol u}_s$. The stress along the surface is related to the surface displacement field by the constitutive relation from Harden & Pleiner (1994), which contains an elastic part proportional to the strain and a viscous part proportional to the strain rate. The surface stress (force per unit length on a differential line element on the surface) $\boldsymbol {\sigma }_s$ is written in terms of the strain and the strain rate fields:

(3.3)\begin{equation} \boldsymbol{\sigma}_s = \left(G_s + \eta_s \frac{\partial}{\partial t} \right) (\nabla_s {\boldsymbol u}_s + (\nabla_s {\boldsymbol u}_s)^{T} - {\frac{1}{2}} {\boldsymbol I} (\nabla_s \boldsymbol{\cdot} {\boldsymbol u}_s)) + \left(B_s + \eta_{bs} \frac{\partial}{\partial t} \right) {\boldsymbol I} (\nabla_s \boldsymbol{\cdot} {\boldsymbol u}_s). \end{equation}

Here, $\nabla _s$ is the gradient operator along the surface of the membrane, $G_s$ and $B_s$ are the surface shear and bulk moduli, which have dimensions of force per unit length, and $\eta _s$ and $\eta _b$ are surface shear and bulk viscosities. At the interface, the normal stress condition is

(3.4)\begin{equation} {\boldsymbol n} \boldsymbol{\cdot} ({\boldsymbol \sigma}_A - {\boldsymbol \sigma}_B) \boldsymbol{\cdot} {\boldsymbol n} = T \nabla_s \boldsymbol{\cdot} {\boldsymbol n}, \end{equation}

where ${\boldsymbol \sigma }_A$ and ${\boldsymbol \sigma }_B$ are the stresses in fluids $A$ and $B$, $T$ is the surface tension, ${\boldsymbol n}$ is the unit normal to the interface and the term on the right-hand side is due to the surface tension of the membrane. The tangential stress balance condition is

(3.5)\begin{equation} ({\boldsymbol I} - {\boldsymbol n} {\boldsymbol n}) \boldsymbol{\cdot} ({\boldsymbol \sigma}_A - {\boldsymbol \sigma}_B) \boldsymbol{\cdot} {\boldsymbol n} = \nabla_s \boldsymbol{\cdot} {\boldsymbol \sigma}_s. \end{equation}

Here, $({\boldsymbol I} - {\boldsymbol n} {\boldsymbol n})$ is the transverse projection operator which projects the stress onto the surface of the membrane. The two stress balance conditions (3.4) and (3.5) are augmented by the continuity of velocity conditions, that is, the velocities of the membrane (partial derivative of displacement with respect to time in the linear approximation) are equal to the fluid velocities. In the viscous limit, the surface displacement model of Thaokar & Kumaran (2002) predicts an instability that is not accessible by the simple membrane model without surface displacement.

3.3. Viscoelastic solid continuum

The standard configurations and coordinate systems used for the flow past viscoelastic continuum surfaces are shown in figure 3. The Couette flow in a channel with one soft wall that is infinite in the direction perpendicular to the plane of the flow (Kumaran, Fredrickson & Pincus 1994; Shankar & Kumaran 2001b; Chokshi & Kumaran 2007, 2008a, 2009) is shown in figure 3(a). Here, the fluid layer of thickness $h_f$ is bounded by a compliant wall of thickness $h_s$ on one side and a rigid wall on the other side, and the rigid wall is moved with a constant velocity to generate flow. A linear velocity profile is generated for a steady fully developed flow. Perturbations are imposed on the interface between the fluid and compliant wall, and the stability of these perturbations is studied. The experimental equivalent of this flow is carried out in a commercial rheometer, as shown in figure 9 and discussed in § 8. A Cartesian coordinate system is used, where $x$ is the flow direction, $y$ is the gradient direction and $z$ is perpendicular to the plane of the flow.

Figure 3. The configurations and representative coordinate systems used for studying the stability of internal flows past soft surfaces. (a) The Couette flow in a channel with one compliant wall, (b) the pressure-driven flow in a channel with two compliant walls and (c) the flow in a compliant tube.

The second configuration is channel flow, where both walls are made of compliant material each of thickness $h_s$, and the fluid has thickness $2 h_f$, as shown in figure 3(b). The flow is driven by a pressure gradient in the streamwise direction. This has been used in Gaurav & Shankar (2010) and Patne & Shankar (2019a,b) for the flow between viscoelastic walls. There are two types of perturbations in this case, the ‘sinuous’ mode where the wall displacement is symmetric about the centreline, and the ‘varicose’ mode where the wall displacement is antisymmetric about the centreline, since an arbitrary perturbation can be expressed as the superposition of the sinuous and varicose modes. The third configuration is the flow in a tube with an annular soft solid forming the wall of the tube shown in figure 3(c) (Kumaran 1995, 1996, 1998; Shankar & Kumar 2004; Shankar & Kumaran 1999, 2000, 2001a, 2002).

In the viscoelastic continuum model, the solid dynamics is characterised by the displacement field, which is the displacement of material points from their reference state due to the stresses exerted on the solid. The reference state has to be selected carefully for the problem under consideration. Consider the steady flow of a fluid past a viscoelastic wall, as shown in figure 4(b). Due to the shear stress exerted by the fluid at the wall, there is elastic deformation of the wall from the equilibrium state with no flow, shown in figure 4(a). Therefore, for the steady deformed state of the solid, the reference state is the equilibrium state with no flow. It is shown in § 3.3.3 that this displacement is unidirectional – the displacement in the streamwise ($x$) direction is only a function of the cross-stream ($y$) direction. The linear stability analysis is carried out for the steady flow and the deformed steady state of the solid, where a perturbation is imposed on the steady state as shown in figure 4(c), and the growth/decay of perturbation is analysed. In this case, it is natural to consider the reference state as the deformed steady state of the solid shown in figure 4(b). In addition to internal displacement of material points within the solid, there could also be a change in the shape of the interface due to the perturbations. The boundary conditions, continuity of velocity and stress, are to be applied at the perturbed interface.

Figure 4. The solid and fluid configurations for (a) the equilibrium state in the absence of flow, (b) steady unidirectional displacement profile (red) in the solid due to a steady flow (blue) in the fluid, (c) perturbed interface in the linear stability analysis and (d) the displacement of material points from the base steady state $(X, Y, Z)$ to the perturbed state $(x,y,z)$ and the unit normal and tangent at the perturbed interface. The red lines in (b) are the steady displacements $\bar {u}_x$ and the red line in (d) is the displacement due to the perturbation $(u_x, u_y, u_z)$.

There are two kinds of descriptions for the solid deformation, the Lagrangian and Eulerian descriptions. Consider a material point at the location ${\boldsymbol X}$ in the reference state which moves to the location ${\boldsymbol x}$ due to deformation of the solid, as shown in figure 5. In the Lagrangian description, the independent coordinate is the location of the material point in the reference state ${\boldsymbol X}$, and the location in the current state is expressed as ${\boldsymbol x}({\boldsymbol X}, t)$. In the Eulerian description, the independent coordinate is the location in the current state, and the location in the reference state prior to application of additional stress is written as ${\boldsymbol X}({\boldsymbol x}, t)$. The additional elastic stress depends on the additional deformation of a material element ${\rm \Delta} {\boldsymbol X}$ in the reference state which becomes ${\rm \Delta} {\boldsymbol x}$ in the current state, as shown in figure 5. The Lagrangian and Eulerian descriptions are described in §§ 3.3.1 and 3.3.2, and the balance laws and constitutive relations are derived. Details of the formulation can be found in the classic texts of Malvern (1969) and Holzapfel (2000). This is followed by the calculation of the deformation at steady state and the linear stability equations using the Lagrangian formulation in §§ 3.3.3 and 3.3.4, respectively, for the flow in a channel. The Lagrangian description is selected because it is easier to implement for reasons discussed at the end of § 3.3.2. The boundary conditions for the continuity of velocity and stress at the interface between the solid and the fluid are derived in § 3.3.5.

Figure 5. The location of material points in the reference state ${\boldsymbol X}$ and ${\boldsymbol X} + {\rm \Delta} {\boldsymbol X}$, in the current state ${\boldsymbol x}$ and ${\boldsymbol x} + {\rm \Delta} {\boldsymbol x}$, and the deformation of a material element ${\rm \Delta} {\boldsymbol X}$ in the reference state to ${\rm \Delta} {\boldsymbol x}$ in the current state.

3.3.1. Lagrangian description

In the Lagrangian description, the coordinate ${\boldsymbol X}$ in the reference configuration is used as the independent coordinate. The material point at the reference location ${{\boldsymbol X} = X {\boldsymbol e}_x + Y {\boldsymbol e}_y + Z {\boldsymbol e}_z}$ moves to the new location ${\boldsymbol x} = x {\boldsymbol e}_x + y {\boldsymbol e}_y + z {\boldsymbol e}_z$ in the deformed state at the time $t$. Here, we follow the common practice of using the same basis vectors for the reference state and the current state. The displacement vector is defined as ${\boldsymbol u} = {\boldsymbol x} - {\boldsymbol X}$:

(3.6)\begin{equation} u_i({\boldsymbol X}, t) = x_i({\boldsymbol X}, t) - X_i.\end{equation}

Here, indicial notation is used to represent vectors, and a repeated index is a dot product. The velocity in the Lagrangian coordinate system is

(3.7)\begin{equation} v_i = \frac{\partial x_i({\boldsymbol X}, t)}{\partial t} = \frac{\partial u_i({\boldsymbol X}, t)}{\partial t}. \end{equation}

The deformation gradient tensor relates the vector distance between two adjacent points in the current state to that in the reference state shown in figure 5, ${\rm \Delta} {\boldsymbol x} = {\boldsymbol F} \boldsymbol {\cdot } {\rm \Delta} {\boldsymbol X}$:

(3.8)\begin{equation} {\rm \Delta} x_i = F_{ij} {\rm \Delta} X_j, \end{equation}

where the deformation gradient tensor is $(\partial {\boldsymbol x}({\boldsymbol X}, t) / \partial {\boldsymbol X})$:

(3.9)\begin{equation} F_{ij} = \frac{\partial x_i}{\partial X_j} = \delta_{ij} + \frac{\partial u_i}{\partial X_j}. \end{equation}

The solid is incompressible if the volume measure of a differential volume element is preserved as the material deforms, or if the Jacobian determinant for the transformation from the reference location ${\boldsymbol X}$ to the deformed location ${\boldsymbol x}$ is $1$. Since the Jacobian matrix is just ${\boldsymbol F}$ (equation (3.8)), the solid is incompressible if

(3.10)\begin{equation} {\textrm{Det}}({\boldsymbol F}) = 1. \end{equation}

The deformation gradient tensor is not invariant under rotation of the form ${\boldsymbol x}' = {\boldsymbol R} \boldsymbol {\cdot } {\boldsymbol x}$. Here, ${\boldsymbol R}$ is the rotation tensor which maps the position ${\boldsymbol x}$ onto a new position ${\boldsymbol x}'$ rotated with respect to ${\boldsymbol x}$; the rotation tensor satisfies the relation ${\boldsymbol R}^{T} = {\boldsymbol R}^{-1}$. A fundamental requirement of any strain measure that is used in the constitutive relation for the stress is that it should be invariant under solid body rotation. This is because solid body rotation does not stretch or compress material line elements, and therefore should not generate a stress. The deformation gradient tensor cannot be used as a strain measure, because it is not invariant under solid body rotation, ${\boldsymbol F} = {\boldsymbol R}$. The left and right Cauchy–Green tensors, ${\boldsymbol F} \boldsymbol {\cdot } {\boldsymbol F}^{T}$ and ${\boldsymbol F}^{T} \boldsymbol {\cdot } {\boldsymbol F}$, respectively, are the simplest strain measures which are invariant under rotation. The left Cauchy–Green tensor is used in the Cauchy and Mooney–Rivlin constitutive relations for the stress defined later in (3.16) and (3.17). In indicial notation, the left Cauchy–Green or finger strain tensor ${\boldsymbol b} = {\boldsymbol F} \boldsymbol {\cdot } {\boldsymbol F}^{T}$ is

(3.11)\begin{equation} b_{ij} = F_{ik} F_{jk} = \delta_{ij} + \frac{\partial u_i}{\partial X_j} + \frac{\partial u_j}{\partial X_i} + \frac{\partial u_i}{\partial X_k} \frac{\partial u_j}{\partial X_k}. \end{equation}

It should be noted that the left Cauchy–Green tensor, (3.11), is a nonlinear function of the gradient of the displacement field, $(\partial u_i/\partial X_j)$. In the ‘linear’ approximation of the deformation gradient tensor, the last term on the right-hand side is neglected in (3.11). When this approximation is made, the strain tensor is not material frame invariant. This approximation can be used only when the measure of the strain is small, $|\partial u_i/\partial X_j| \ll 1$. The linear approximation is valid for the progression from the steady state (figure 4b) to the perturbed state (figure 4c), since the perturbation amplitudes are considered infinitesimal. This approximation is not valid, in general, for the progression from the equilibrium state (figure 4a) to the steady state (figure 4b). It is argued at the end of § 3.3.5 that the linearisation approximation is not valid for pressure-driven flow at low Reynolds number, and the neglect of the nonlinear term in (3.11) can qualitatively change the result of the stability analysis. The linearisation approximation is valid at high Reynolds number, where the strain in the base state is small.

The time derivative $\dot {{\boldsymbol F}}$ of ${\boldsymbol F}$ is

(3.12)\begin{equation} \dot{F_{ij}} = \frac{\partial}{\partial t} \left( \frac{\partial x_i}{ \partial X_j} \right) = \frac{\partial v_i}{\partial X_j}. \end{equation}

Note that the gradient here is with respect to the reference coordinates and so it is not the spatial velocity gradient. To express this in terms of the spatial velocity gradient, the time derivative $\dot {{\boldsymbol F}}$ can be written as

(3.13)\begin{equation} \dot{F_{ij}} = \frac{\partial x_k}{\partial X_j} \frac{\partial v_i}{\partial x_k} = \frac{\partial v_i}{\partial x_k} F_{kj} = l_{ik} F_{kj}, \end{equation}

where ${\boldsymbol v} = (\partial {\boldsymbol x}({\boldsymbol X},t)/\partial t)$ is the local velocity. The spatial velocity gradient is defined as $l_{ij} = (\partial v_i / \partial x_j)$; this is the rate of deformation tensor in fluid mechanics. The rate of deformation tensor can be written in terms of the deformation gradient tensor as ${\boldsymbol l} = \dot {{\boldsymbol F}} \boldsymbol {\cdot } {\boldsymbol F}^{-1}$:

(3.14)\begin{equation} l_{ij} = \dot{F}_{ik} F_{kj}^{{-}1}. \end{equation}

The constitutive relation for the Cauchy stress tensor typically contains the pressure to satisfy incompressibility, the elastic stress due to the deformation gradient and the viscous stress due to the strain rate. In its most general form, the elastic stress is a function of the left Cauchy–Green tensor ${\boldsymbol b}$ (equation (3.11)), which is a symmetric tensor. In the theory of hyperelastic materials, the stress is written as the functional derivative of a free energy $\psi (I_1, I_2, I_3)$, where $I_1 = \textrm {Tr}{{\boldsymbol b}}$, $I_2 = \textrm {Tr}{{\boldsymbol b}^{-1}} \textrm {Det}{{\boldsymbol b}}$ and $I_3 = \textrm {Det}{{\boldsymbol b}}$ are the three principal scalar invariants of ${\boldsymbol b}$. Since $\textrm {Det}{{\boldsymbol b}} = 1$ is a constant for an incompressible elastic material (see (3.10)), the relevant invariants are $I_1$ and $I_2$. The elastic stress is the derivative of the free energy with respect to the right Cauchy–Green tensor ${\boldsymbol c} = {\boldsymbol F}^{T} \boldsymbol {\cdot } {\boldsymbol F}$, $\boldsymbol {\sigma } = (\partial \psi /\partial {\boldsymbol c})$, which is written in indicial notation as

(3.15)\begin{equation} \sigma_{ij}^{e} = 2 \left( \frac{\partial \psi}{\partial I_1} b_{ij} - \frac{\partial \psi}{\partial I_2} b^{{-}1}_{ij} \right), \end{equation}

where $b^{-1}_{ij} = ({\boldsymbol b}^{-1})_{ij}$ is the $i,j$ element of the inverse of matrix ${\boldsymbol b}$. The details of this derivation are provided in Malvern (1969) and Holzapfel (2000). The neo-Hookean model for the elastic stress is obtained by substituting $\psi = G I_1$ in (3.15):

(3.16)\begin{equation} \sigma_{ij}^{e} = G b_{ij}, \end{equation}

where $G$ is the elasticity modulus. The more general Mooney–Rivlin model is obtained by substituting $\psi = \frac {1}{2} (G_1 I_1 - G_2 I_2)$ in (3.15):

(3.17)\begin{equation} \sigma_{ij}^{e} = G_1 b_{ij} + G_2 b^{{-}1}_{ij}. \end{equation}

More complicated forms of the free energy can be used to derive equations for the stress that are functions of ${\boldsymbol b}$ and ${\boldsymbol b}^{-1}$, and the two invariants of ${\boldsymbol b}$. However, (3.17) is the simplest relationship in which the stress tensor is a linear function of the strain tensor ${\boldsymbol b}$ or its inverse.

The constitutive relation for the stress tensor also contains an isotropic pressure required to satisfy incompressibility, and a viscous part due to the velocity gradient, (3.14). The latter is considered proportional to the symmetric part of the velocity gradient, analogous to Newton's law of viscosity for Newtonian fluids. The constitutive relation for the stress tensor is

(3.18)\begin{equation} \sigma^{s}_{ij} ={-} p_s \delta_{ij} + G b_{ij} + \mu_s (l_{ij} + l_{ji}),\end{equation}

where $p_s$ is the dynamical pressure in the solid and $\mu _s$ is the solid viscosity. The first term on the right-hand side in (3.18) is the isotropic pressure required to satisfy incompressibility, the second term is the elastic stress due to the deformation gradient and the third term is the viscous stress proportional to the symmetric part of the velocity gradient tensor.

The momentum conservation equation in the Lagrangian reference frame is

(3.19)\begin{equation} \rho_s \frac{\textrm{D} v_i}{\textrm{D} t} = \frac{\partial P_{ij}}{\partial X_j}, \end{equation}

where $(\textrm {D} / \textrm {D} t)$ is the substantial derivative $(\partial / \partial t + {\boldsymbol v} {\boldsymbol \cdot } \boldsymbol {\nabla })$. The Piola–Kirchoff stress tensor is ${\boldsymbol P} = ({\boldsymbol F}^{-1} \boldsymbol {\cdot } {\boldsymbol \sigma }_s)^{T}$,

(3.20)\begin{equation} P_{ij} = F^{{-}1}_{jk} \sigma^{s}_{ki}, \end{equation}

which is the stress tensor when expressed in terms of the area element and unit normal in the reference coordinate system. It is necessary to pre-multiply the stress tensor ${\boldsymbol \sigma }_s$ by ${\boldsymbol F}^{-1}$ in the Lagrangian description, because all the derivatives are with respect to the reference coordinates. The Piola–Kirchoff stress tensor is obtained by projecting the local surface force (stress dotted with unit normal) at the deformed location onto the original location. It should be noted that the stress tensor ${\boldsymbol \sigma }_s$ in the spatial coordinates is symmetric, and it is identical to the stress tensor used in fluid dynamics. The Piola–Kirchoff stress tensor in the reference coordinates is not symmetric in general. For the neo-Hookean constitutive relation (3.18), the momentum conservation equation for the Lagrangian formulation is

(3.21)\begin{equation} \rho_f \frac{\textrm{D} v_i}{\textrm{D} t} ={-} \frac{\partial F^{{-}1}_{ji} p_s}{\partial X_j} + G \frac{\partial F_{ij}}{\partial X_j} + \mu_s \left( \frac{\partial (F^{{-}1}_{jk} (l_{ik} + l_{ki}))}{\partial X_j} \right). \end{equation}

The reference state has to be carefully defined in the Lagrangian description. The definition of the reference state is not important when the linear model is used for the relationship between the stress and displacement gradient, but it is of critical importance when ${\boldsymbol b}$ (equation (3.11)) is a nonlinear function of the gradients in the displacement field. In a channel/tube flow with compliant walls, there is a fluid stress exerted at the interface between the fluid and the wall, which deforms the solid. For the base state in which there is a steady flow (figure 4b), the reference state is the equilibrium (zero-stress) state in the absence of flow (figure 4a). The solid stress in the base state is calculated using a Lagrangian description in which the position coordinates are those in the reference state. This is illustrated in § 3.3.3. For the linear stability calculations in the perturbed state (figure 4c), the state of the solid deformed due to the steady flow (figure 4b) is chosen as the reference state, as shown in § 3.3.4.

The above formulation has been used in Chokshi & Kumaran (2007, 2008a,b), and is referred to as the L3 formulation in Patne, Giribabu & Shankar (2017) and Patne & Shankar (2019a). This is because three states, the equilibrium (no-flow) state, the stressed steady state and the perturbed state, are considered in the analysis. An alternative formulation has been used by Gkanis & Kumar (2005), Gaurav (2009) and Gaurav & Shankar (2010) where there are only two states; the reference state for the linear stability analysis is considered as the unstressed equilibrium state for the solid. This is referred to as the L2 model in Patne et al. (2017) and Patne & Shankar (2019a), where it is shown that the errors in the Taylor series expansion could be significant. This is because the solid displacement from the reference equilibrium state to the stressed base state may not be small. Thus, it is best to use the formulation in §§ 3.3.3 and 3.3.4, denoted the L3 formulation in Patne et al. (2017) and Patne & Shankar (2019a).

The interfacial boundary conditions for the velocity and stress in the Lagrangian description are applied as follows. On the solid side, the values of the velocity and stress at the reference location are used in the analysis. This is because the position label for the particle velocity and stress in the Lagrangian description is still the reference location, and not the current location. Thus, there is no necessity to use a Taylor expansion of the quantities at the perturbed location about their values at the unperturbed location. In contrast, the fluid is described in an Eulerian reference frame. In this case, it is necessary to use a Taylor expansion of the quantities at the perturbed interface location about those in the base state.

3.3.2. Eulerian description

In the Eulerian description, the displacement field is expressed in terms of the current coordinate, ${\boldsymbol x}$. If we consider a material point with unstressed coordinate ${\boldsymbol X}$ which moves to the location ${\boldsymbol x}$ upon application of a stress, the new location is expressed as

(3.22)\begin{equation} X_i({\boldsymbol x}, t) = x_i - u_i({\boldsymbol x}, t). \end{equation}

where ${\boldsymbol u}$, the displacement field, is a function of the current position ${\boldsymbol x}$. In this case, the relation between a differential displacement in the deformed and reference coordinates is

(3.23)\begin{equation} f_{ij} = \frac{\partial X_i}{\partial x_j} = \delta_{ij} - \frac{\partial u_i}{\partial x_j}, \end{equation}

where $\delta _{ij}$ is the identity tensor. Comparing (3.9) and (3.23), it is evident that ${\boldsymbol f} = {\boldsymbol F}^{-1}$. Therefore, mass conservation equation is

(3.24)\begin{equation} {\textrm{Det}}({\boldsymbol f}) = 1 \end{equation}

or

(3.25)\begin{equation} \frac{\partial v_i}{\partial x_i} = 0, \end{equation}

where the velocity field ${\boldsymbol v}$ is given by

(3.26)\begin{equation} v_i = \left( \frac{\partial x_i}{\partial t} \right)_{{\boldsymbol X}} = \frac{\partial u_i}{\partial t} + v_j \frac{\partial u_i}{\partial x_j}.\end{equation}

Thus, the velocity field in the Eulerian description has to be obtained from an implicit relation between the time derivative of the displacement field and its spatial gradients.

The elastic stresses for the neo-Hookean and Mooney–Rivlin models are given by (3.16) and (3.17), with ${\boldsymbol F} = {\boldsymbol f}^{-1}$. Here, the constitutive relations are equivalent only if we make the substitution ${\boldsymbol F} \boldsymbol {\cdot } {\boldsymbol F}^{T} \rightarrow ({\boldsymbol f}^{T} \boldsymbol {\cdot } {\boldsymbol f})^{-1}$ and $({\boldsymbol F} \boldsymbol {\cdot } {\boldsymbol F}^{T})^{-1} \rightarrow {\boldsymbol f}^{T} \boldsymbol {\cdot } {\boldsymbol f}$. Similarly, the rate of deformation tensor in the Eulerian framework is $l_{ij} = (\partial v_i/\partial x_j)$, where the velocity ${\boldsymbol v}$ is determined from (3.26). The stress tensor for the neo-Hookean model in the Eulerian formulation is

(3.27)\begin{equation} \sigma^{s}_{ij} ={-} p_s \delta_{ij} + G ({\boldsymbol f}^{T} \boldsymbol{\cdot} {\boldsymbol f})^{{-}1} + \mu_s (l_{ij} + l_{ji}). \end{equation}

The mass conservation equation is given by (3.24) and the momentum conservation equation is

(3.28)\begin{equation} \rho \frac{\textrm{D} v_i}{\textrm{D} t} = \frac{\partial \sigma^{s}_{ij}}{\partial x_j}. \end{equation}

In the Eulerian description, the interface conditions are enforced at the perturbed surface on the fluid and solid side. This is accomplished using a Taylor expansion of the quantities at the perturbed surface about their values at the unperturbed surface, in a manner similar to that for the fluid. For a linear stability analysis, only the terms linear in the perturbation quantities are retained, but higher-order terms are also included in weakly nonlinear studies.

There are advantages and disadvantages to both descriptions of solid deformation. The Eulerian description is similar to the flow equations, and might be considered the preferred candidate for the fluid–solid interaction. The disadvantages are that the velocity field (and therefore the rate of deformation tensor) in the Eulerian description is determined from an implicit relation, (3.26), and that the boundary conditions are imposed at the perturbed location of the interface. Another advantage of the Eulerian formulation is that the linearisation of the continuity conditions at the interface is unambiguous, and this can be used to check the consistency of the Taylor expansion of the interface conditions in the Lagrangian formulation. In contrast, the velocity and rate of deformation tensor for the Lagrangian description are determined from explicit relations, (3.7) and (3.14). In addition, the boundary conditions are imposed at the locations in the reference state, and not the current state, since the positions in the Lagrangian description are labelled by their coordinates in the reference state. For these reasons, the linear stability analysis is easier to implement in the Lagrangian formulation compared with the Eulerian formulation. The steady deformation for the flow in a channel is calculated in § 3.3.3 and the equations for the linear stability analysis are derived in § 3.3.4.

3.3.3. Solid deformation at steady state

The configuration, shown in figure 3(b), consists of a layer of fluid of thickness $h_f$ in the region $0 < y < h_f$ adjacent to a viscoelastic solid of thickness $h_s$ in the region $-h_s < y < 0$. In the base state, the fluid flow is a unidirectional flow in the $x$ direction. A general velocity profile with pressure gradient $(dp_f/{\textrm {d} x})$, zero velocity and strain rate $\dot {\gamma }_w$ at the wall is considered. It should be noted that the solid is stressed but stationary in the base state. The deformation gradient tensor ${\boldsymbol F}$ and its inverse are

(3.29a,b) \begin{equation} {\boldsymbol F} = \begin{pmatrix} 1 + \dfrac{\partial u_x}{\partial X} & \dfrac{\partial u_x}{\partial Y} \\ \dfrac{\partial u_y}{\partial X} & 1 + \dfrac{\partial u_y}{\partial Y} \end{pmatrix},\quad {\boldsymbol F}^{{-}1} = \begin{pmatrix} 1 + \dfrac{\partial u_y}{\partial Y} & - \dfrac{\partial u_x}{\partial Y} \\ - \dfrac{\partial u_y}{\partial X} & 1 + \dfrac{\partial u_x}{\partial X} \end{pmatrix}. \end{equation}

In the Lagrangian description, the components of the momentum conservation equation for the steady displacement field $(\bar {u}_x, \bar {u}_y)$ (equation (3.21)), are

(3.30)$$\begin{gather} - \frac{\partial (\bar{p}_s \bar{F}_{xx}^{{-}1})}{\partial X} - \frac{\partial (\bar{p}_s \bar{F}_{yx}^{{-}1})}{\partial Y} + G \left( \frac{\partial \bar{F}_{xx}}{\partial X} + \frac{\partial \bar{F}_{xy}}{\partial Y} \right) = 0, \end{gather}$$

(3.31)$$\begin{gather}- \frac{\partial (\bar{p}_s \bar{F}_{xy}^{{-}1})}{\partial X} - \frac{\partial (\bar{p}_s \bar{F}_{yy}^{{-}1})}{\partial Y} + G \left( \frac{\partial \bar{F}_{yx}}{\partial X} + \frac{\partial \bar{F}_{yy}}{\partial Y} \right) = 0. \end{gather}$$

Here, the overbars denote quantities in stressed base state. The above equations are satisfied when the displacement field is unidirectional along the flow direction, that is, $\bar {u}_x = \bar {u}_x(Y)$ is independent of $X$ and $\bar {u}_y = 0$. In this case, the stress tensor is

(3.32) \begin{equation} \bar{\boldsymbol{\sigma}}_s = \begin{pmatrix} - \bar{p}_s + G \left( 1 + \left(\dfrac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \right)^{2} \right) & G \dfrac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \\ G \dfrac{\textrm{d} \bar{u}_x}{\textrm{d} Y} & - \bar{p}_s + G \end{pmatrix}.\end{equation}

The components of the simplified momentum equation, (3.30) and (3.31), are

(3.33)$$\begin{gather} {-}\frac{\partial \bar{p}_s}{\partial X} + G \frac{\textrm{d}^{2} \bar{u}_x}{\textrm{d} Y^{2}} = 0, \end{gather}$$

(3.34)$$\begin{gather} {-}\frac{\partial \bar{p}_s}{\partial Y} + \frac{\partial \bar{p}_s}{\partial X} \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} = 0. \end{gather}$$

The above equations are solved subject to the zero-displacement boundary conditions $\bar {u}_X = 0$ at $Y=-h_s$, and the stress balance conditions tangential and normal to the surface at the interface $Y=0$:

(3.35a,b)\begin{equation} G \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} = \mu_f \dot{\gamma}_w \quad \text{and} \quad \bar{p}_s + G = \bar{p}_f, \end{equation}

respectively. Here, $\dot {\gamma }_w$ is the fluid strain rate at the wall and $\bar {p}_f$ is the mean pressure in the fluid, which is a linear function of $X$. The mean displacement and pressure fields in the solid are

(3.36)$$\begin{gather} \bar{u}_x ={-} \frac{h_s^{2}}{2 G} \frac{\textrm{d} \bar{p}}{\textrm{d} x} \left( 1 - \frac{Y^{2}}{h_s^{2}} \right) + \frac{\mu_f \dot{\gamma}_w (Y+h_s)}{G}, \end{gather}$$

(3.37)$$\begin{gather}\bar{p}_s = \bar{p} + \left(\frac{\textrm{d} \bar{p}}{\textrm{d} x} \right)^{2} \frac{Y^{2}}{2 G} + \frac{\textrm{d} \bar{p}}{\textrm{d} x} \frac{\mu_f \dot{\gamma}_w Y}{G} - G. \end{gather}$$

The steady displacement field in a tube can be derived in cylindrical coordinates by a procedure similar to that used here (Gaurav 2009) for the neo-Hookean wall model. Here, the displacement field in the solid is unidirectional along the axis, and the axial displacement field is only a function of the radial coordinate.

3.3.4. Linear stability equations

The displacement and pressure fields in the solid are expressed in a manner similar to the velocity field in (2.6):

(3.38a,b)\begin{equation} u_i = \bar{u}_i + \tilde{u}_i \exp{(\imath k_x X + \imath k_z Z + s t)}, \quad p_s = \bar{p}_s + \tilde{p}_s \exp{(\imath k_x X + \imath k_z Z + s t)}, \end{equation}

where $s$ is the growth rate and $k_x, k_z$ are the wavenumbers in the $x$ and $z$ directions. When these are substituted into the mass and momentum equations, (3.10) and (3.21), and linearised in the perturbations $\tilde {u}_i$ and $\tilde {p}_s$, we obtain

(3.39)\begin{align} &\hspace{6.5pc}\frac{\partial \tilde{u}_y}{\partial Y} + \imath k_x \tilde{u}_x + \imath k_z \tilde{u}_z - \imath k_x \tilde{u}_y \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} = 0, \end{align}

(3.40)\begin{align} \rho_s s^{2} \tilde{u}_x &={-} \imath k_x \tilde{p}_s - \frac{\partial \bar{p}_s}{\partial X} \left( \frac{\textrm{d} \tilde{u}_y}{\textrm{d} Y} + \imath k_z \tilde{u}_z \right) + \frac{\partial \bar{p}_s}{\partial Y} \imath k_x \tilde{u}_y + (G + \mu_s s) \left( \frac{\textrm{d}^{2}}{\textrm{d} Y^{2}} - k_x^{2} - k_z^{2} \right) \tilde{u}_x \nonumber\\ & \quad - \mu_s s \left( \tilde{u}_x \left( \imath k_x \frac{\textrm{d}^{2} \bar{u}_x}{d Y^{2}} + k_x^{2} \left( \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \right)^{2} \right) + 2 \imath k_x \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \frac{\textrm{d} \tilde{u}_x}{\textrm{d} Y} \right), \end{align}

(3.41)\begin{align} \rho_s s^{2} \tilde{u}_y &={-} \frac{\textrm{d} \tilde{p}_s}{\textrm{d} Y} + \imath k_x \tilde{p}_s \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} - \imath (k_x \tilde{u}_x + k_z \tilde{u}_z) \frac{\partial \bar{p}_s}{\partial Y} + \frac{\partial \bar{p}_s}{\partial X} \left( \frac{\textrm{d} \tilde{u}_x}{\textrm{d} Y} + \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \imath k_z \tilde{u}_z \right) \nonumber\\ & \quad + (G + \mu_s s) \left( \frac{\textrm{d}^{2}}{\textrm{d} Y^{2}} - k_x^{2} - k_z^{2} \right) \tilde{u}_y \nonumber\\ & \quad - \mu_s s \left(\tilde{u}_y \left(\imath k_x \frac{\textrm{d}^{2} \bar{u}_x}{d Y^{2}} + k_x^{2} \left( \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \right)^{2} \right) + 2 \imath k_x \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \frac{\textrm{d} \tilde{u}_y}{\textrm{d} Y} \right), \end{align}

(3.42)\begin{align} \rho_s s^{2} \tilde{u}_z &={-} \imath k_z \tilde{p}_s + \frac{\partial \bar{p}_s}{\partial X} \imath k_z \left(\tilde{u}_x - \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \tilde{u}_y \right) + \frac{\partial \bar{p}_s}{\partial Y} \imath k_z \tilde{u}_y + (G + \mu_s s) \left( \frac{\textrm{d}^{2}}{\textrm{d} Y^{2}} - k_x^{2} - k_z^{2} \right) \tilde{u}_z \nonumber\\ & \quad - \mu_s s \left( \tilde{u}_z \left(\imath k_x \frac{\textrm{d}^{2} \bar{u}_x}{\textrm{d} Y^{2}} + k_x^{2} \left( \frac{\textrm{d} \bar{u}_X}{\textrm{d} Y} \right)^{2} \right) + 2 \imath k_x \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \frac{\textrm{d} \tilde{u}_z}{\textrm{d} Y} \right). \end{align}

It is evident that though (3.39)–(3.42) for the solid displacement field have a structure similar to (2.7)–(2.10) for the fluid velocity field, the former are more complex. This is due to the nonlinear dependence of the strain tensor on the displacement field in a rotational frame invariant formulation.

3.3.5. Interface conditions

The boundary conditions at the interface between the fluid and the solid are the continuity of velocity and the continuity of tangential and normal stress. The fluid velocity and stress are calculated at the perturbed interface location. On the solid side, the procedure for applying the boundary conditions depends on the formulation. As stated above, for the Lagrangian formulation, the boundary condition is applied at the location $(X, Y=0, Z)$ in the base state.

In the linear approximation, the velocity and the stress in the fluid at the deformed interface location can be separated into two parts, the first due to the perturbations and the second due to the difference in the mean quantities between the perturbed and base interface location. After carrying out a Taylor expansion for the fluid velocity about the base state, the following velocity boundary conditions are applicable at the location $y=Y=0$ of the interface in the base state:

(3.43a,b)\begin{equation} s \tilde{u}_x = \tilde{v}_x + \tilde{u}_y \dot{\gamma}_w \quad \text{and} \quad s \tilde{u}_y = \tilde{v}_y. \end{equation}

In the tangential velocity balance condition, the fluid velocity is the sum of the velocity perturbation $\tilde {v}_x$ and the mean fluid velocity at the deformed interface which is proportional to the perturbation $\tilde {u}_y$ times the fluid strain rate $\dot {\gamma }_w$ at the undisturbed surface. There is no contribution analogous to the latter on the solid side, because the solid strain rate is zero in the base state. The perturbations to the fluid velocity and solid displacement in (3.43a,b) are calculated at the undeformed interface in the linear approximation, since the correction due to deformation is quadratic in the perturbation amplitudes.

The normal stress balance condition at the location $y=Y=0$ is

(3.44)\begin{equation} - \tilde{p}_s + 2 \left(G + \mu_s s \right) \frac{\textrm{d} \tilde{u}_y}{\textrm{d} Y} - 2 \mu_s s \imath k_x \tilde{u}_y \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} ={-} \tilde{p}_f + 2 \mu_f \frac{\textrm{d} \tilde{v}_y}{\textrm{d} y} + \tilde{u}_x \frac{\textrm{d} \bar{p}}{\textrm{d} x}. \end{equation}

In the above equation, there is no variation in the normal stress in the fluid in the base state, and so there is no equivalent of the second term on the right-hand side in the tangential velocity balance condition in (3.43a,b). However, there is a variation in the normal stress in the fluid in the streamwise direction due to the mean pressure gradient, resulting in the last term on the right-hand side in the normal stress equation. Though there is a variation in the normal stress in the solid, there is no contribution in (3.44) due to this variation, because the normal stress boundary condition is applied at the base state location in the Lagrangian formulation.

The balance equation for the tangential stress in the streamwise direction along the surface is

(3.45)\begin{align} & ( G + \mu_s s ) \left( \frac{\textrm{d} \tilde{u}_x}{\textrm{d} Y} + \imath k_x \tilde{u}_y \right) + G \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \frac{\textrm{d} \tilde{u}_y}{\textrm{d} Y} - \mu_s s \imath k_x \tilde{u}_x \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} - G \tilde{u}_y \left( \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} \right)^{2} \nonumber\\ &\quad = \mu_f \left( \frac{\textrm{d} \tilde{v}_x}{\textrm{d} y} + \frac{\textrm{d} \tilde{v}_y}{\textrm{d} x} \right) + \mu_f \tilde{u}_y \frac{\textrm{d}^{2} \bar{v}_x}{{\textrm{d} y}^{2}}. \end{align}

The first term on the right-hand side in (3.45) is the fluid shear stress due to the velocity perturbation evaluated at the unperturbed surface. The second term on the right-hand side is the product of the surface displacement and the gradient in the mean fluid shear stress, $\mu _f (\textrm {d} \bar {v}_x/{\textrm {d} y})$, at the surface. The first three terms on the left-hand side are the elastic and viscous shear stress perturbation at the reference location. The last term on the left-hand side in (3.45) is due to the variation to the tangent in the surface caused by the surface displacement. The unit normal ${\boldsymbol n}$ and tangent ${\boldsymbol t}$ in the perturbed state for the channel geometry (figure 4d) are

(3.46a,b)\begin{equation} {\boldsymbol n} = {\boldsymbol e}_y - \frac{\textrm{d} u_y}{\textrm{d} x} {\boldsymbol e}_x, \quad {\boldsymbol t} = {\boldsymbol e}_x + \frac{\textrm{d} u_y}{\textrm{d} x} {\boldsymbol e}_y, \end{equation}

where $u_y$ is the normal surface displacement. The expressions in (3.46a,b) have been linearised in the perturbation amplitude $u_y$. The stress tangential to the surface on the solid side, ${\boldsymbol t} \boldsymbol {\cdot } {\boldsymbol \sigma }_s \boldsymbol {\cdot } {\boldsymbol n}$, is $\tilde {\sigma }_{xy}^{s} + ({\textrm {d} \tilde {u}_y}/{\textrm {d} x}) (\bar {\sigma }^{s}_{yy} - \bar {\sigma }^{s}_{xx})$ in the linear approximation for the surface displacement. From (3.32), the normal stress difference $\bar {\sigma }_{xx} - \bar {\sigma }_{yy} = G (\textrm {d} \bar {u}_x/ \textrm {d} Y)^{2}$; this leads to the last term on the left-hand side in (3.46a,b). Since there is no normal stress difference in the base state for a Newtonian fluid, there is no equivalent on the right-hand side.

Finally, the balance equation for the tangential stress perpendicular to the plane of flow is

(3.47)\begin{equation} (G + \mu_s s) \left(\frac{\textrm{d} \tilde{u}_z}{\textrm{d} Y} + \imath k_z \tilde{u}_y \right) - \mu_s s \imath k_x \tilde{u}_z \frac{\textrm{d} \bar{u}_x}{\textrm{d} Y} = \mu_f \left(\frac{\textrm{d} \tilde{v}_z}{\textrm{d} y} + \imath k_z \tilde{v}_y \right). \end{equation}

The magnitude of the solid strain in the base state can be estimated based on the fluid stress acting at the interface. The fluid stress scales as $(\mu _f V_f / h_f)$, and therefore the solid strain in the base state is $(\mu _f V_f /G h_f)$, which is the dimensionless parameter $\varGamma$. For a viscous flow in the limit of low Reynolds number (${Re} \ll 1$), the analysis of § 5 indicates that there is an instability when the parameter $\varGamma$ exceeds a critical value. In this case, the solid strain is $O(1)$, and it is necessary to retain the mean strain in the perturbation (3.39)–(3.42) for the solid displacement field, and the boundary conditions (3.43a,b)–(3.45).

3.3.6. Incompressibility

The stability studies have all assumed that the solid is incompressible ((3.10) and (3.24)) and the solid displacement field in the base state is unidirectional. This facilitates an expansion of the velocity field in Fourier modes in the streamwise direction. In real applications, materials such as polymer gels are compressible, but their compression modulus is $10 - 10^{2}$ times larger than the shear modulus (Verma & Kumaran 2012; Srinivas & Kumaran 2017b). The incompressibility approximation is a good one for flows driven by moving boundaries, such as a Couette flow, but may not be valid for pressure-driven flows. In a pressure-driven flow, the pressure difference across the length of the tube/channel is larger than the wall shear stress by a factor $(L/h_f)$, where $L$ is the length of the tube/channel and $h_f$ is the characteristic cross-stream dimension. For sufficiently long conduits, the pressure difference could be significantly larger than the wall shear stress, and the compression could be comparable to or larger than the shear strain. The effect of wall deformation has not been considered so far in theoretical studies, though experimental studies (Verma & Kumaran 2012, 2013; Srinivas & Kumaran 2017b) show that this could be an important effect. A ‘local’ stability analysis was performed by Verma & Kumaran (2013), where the local conduit dimension and pressure gradient at different downstream distances were used in a linear analysis to determine the stability of the system. The implicit assumption here is that the wavelength of the perturbations is small compared to the characteristic length for flow evolution. With this assumption, the results of the linear stability analysis are found to be in reasonable agreement with experimental results.

4. Non-dimensionalisation and numerical methods

The dimensional parameters in the problem are as follows. The length scales are the thickness of the fluid $h_f$ and the solid $h_s$, and the velocity scale is the characteristic fluid velocity $V_f$. The fluid properties are the density $\rho _f$ and viscosity $\mu _f$, and the solid properties are the density $\rho _s$, the modulus of elasticity $G$ and the viscosity $\mu _s$. In the linear stability analysis, the parameters are the wavenumber $k$ and the growth rate $s$. In addition to the Reynolds number ${Re} = (\rho _f V_f h_f /\mu _f)$, there is one other dimensionless group involving the solid elasticity modulus. This is usually expressed as the velocity-independent parameter $\varSigma = (\rho _f G h_f^{2}/\mu _f^{2})$, though sometimes the ratio of viscous and elastic stress, $\varGamma = (V_f \mu _f/G h_f)$, is also used. The other dimensionless parameters are the ratio of the fluid and solid viscosities $\mu _r = \mu _s/\mu _f$ and the ratio of the characteristic lengths of the fluid and solid $H = (h_s / h_f)$. In § 5, the limit of zero Reynolds number is considered, where there is a balance between the viscous stress in the fluid and the elastic stress in the solid for $\varGamma \sim 1$. In this case, the flow stability depends only on the dimensionless parameters $\varGamma$, $\mu _r$ and $H$. At high Reynolds number, the critical Reynolds number for the onset of instability scales as a power of the parameter $\varSigma$. It is shown in § 6.1 that ${Re}_c \propto \varSigma ^{1/2}$ for the high-Reynolds-number modes that have an internal critical layer within the flow, and in § 6.3 it is shown that ${Re} \propto \varSigma ^{3/4}$ for the wall modes.

It is possible to adopt a uniform non-dimensionalisation for some parameters. The spatial coordinates are all non-dimensionalised by the characteristic length $h_f$ and the wavenumber by $h_f^{-1}$. For the solid dynamics, the pressure and the elastic stress are non-dimensionalised by the elasticity modulus $G$. The superscript $^{\ast }$ is used for quantities non-dimensionalised in this way. The non-dimensionalisation of the fluid velocity and pressure, the growth rate and the wave velocity are different in the low- and high-Reynolds-number regimes. At low Reynolds number, the time dimension is non-dimensionalised by the characteristic time scale $(\mu _f/G)$; quantities non-dimensionalised in this manner are denoted by the superscript $^{\ast }$. At high Reynolds number, $(h_f/V_f)$ is used to non-dimensionalise the time dimension, and the superscript $^{{{\dagger}} }$ is used for quantities non-dimensionalised in this manner. The scalings for the different quantities are summarised in table 1.

Table 1. (a) The non-dimensionalisation used for the low-Reynolds-number analysis in § 5 and high-Reynolds-number analysis in § 6 for the flow past continuum viscoelastic surfaces and the scalings for the properties of the spring-backed wall in § 6.1. (b) The dimensionless groups for the continuum viscoelastic wall model.

The linear stability equations reduce to ordinary differential equations in the cross-stream coordinate, (2.7)–(2.10) for the fluid and (3.39)–(3.42) for the wall made of a viscoelastic continuum. The solid and fluid equations are solved separately in their separate domains. The boundary conditions at the perturbed interface are expressed using a Taylor series expansion about the unperturbed interface in (3.43a,b)–(3.47). These are solved either using numerical integration in the cross-stream coordinate or using spectral methods.

In the spatial integration scheme, the equations for the solid and fluid are numerically integrated in their respective domains. Linearly independent solutions that satisfy the boundary conditions at the fixed boundaries are determined. The fixed boundaries are, for example, $y = h_f$ for the fluid and $Y = - h_s$ for the solid for the Couette flow in figure 3(a), while symmetry conditions are applied, for example, at the centre of the tube $r = 0$ in figure 3(c). At the interface, the boundary conditions are expressed as the product of a coefficient matrix multiplying the column vector of the coefficients of the linearly independent solutions. The growth rate is determined from the condition that the determinant of the coefficient matrix is zero. Note that the elements of the coefficient matrix are complex in general, so the real and imaginary parts of the growth rate are determined from the zero determinant condition. Typically, the simulation starts with an initial guess for the growth rate, and a Newton–Raphson iteration scheme is used to converge to the solution.

The second method is the spectral method (Orszag 1971; Boyd 2000), where the displacement and velocity fields are usually expressed as the sum of a finite set of $N$ Chebyshev polynomials each. These are substituted into the four governing equations, two each for the solid and fluid, and which are set equal to zero at $(4N - 8)$ Gauss–Labatto collocation points in the solid and fluid domains. The boundary conditions, two each at $y = h_f$ and $Y = - h_s$, and the four interface conditions provide an additional eight rows in the pseudospectral matrix. This results in a polynomial eigenvalue problem of dimension $4N \times 4N$ for the growth rate $s$, which is solved numerically. The number of solutions for the growth rate is equal to the dimension of the matrix. The physical solutions are identified as those which converge as the number of basis functions $N$ is increased.

The spatial integration method is more accurate, since a Newton–Raphson method is used to obtain quadratic convergence. However, this method requires an initial guess for the eigenvalue, and is difficult to use when there are multiple solutions for the eigenvalue. The spectral method has relatively poorer accuracy, since the solution depends on the number of basis functions used. It is preferable to use a composite scheme, where the eigenvalues are first estimated using the spectral scheme, and these are used as the initial guess in a spatial integration scheme to verify that they are not spurious, and to improve the accuracy.

5. Viscous flow

The stability analysis is considerably simplified for zero Reynolds number, where the inertial terms in (2.8)–(2.10) for the fluid and in (3.40)–(3.42) for the solid are neglected. Though the fluid equations are quasi-steady, there is time dependence in the equations for the solid for $\mu _r \neq 0$, because the velocity field is the time derivative of the displacement field in the solid equations. When the solid viscosity is zero, the equations for the solid are also quasi-steady and have no explicit time dependence. However, the time dependence enters through the interface velocity continuity, where the time derivative of the displacement in the solid is equal to the velocity in the fluid. This interface coupling can cause a flow instability even in the absence of inertia for a continuum viscoelastic wall model. The zero-Reynolds-number analysis for the Couette flow is first discussed, and a sample calculation is presented to illustrate the simplifications in the viscous limit. The stability of a pressure-driven flow is then discussed. For a viscous flow, it is appropriate to scale time by $(\mu _f/G)$ and velocity by $(h_f G / \mu _f)$. The scaled velocity is defined as $\tilde {v}^{\ast }_i = (\mu _f \tilde {v}_i / h_f G)$ and the scaled growth rate is $s^{\ast } = (s \mu _f/G)$. The non-dimensional fluid pressure is $p_f^{\ast } = (p_f/G)$. The other quantities are non-dimensionalised as shown in table~1.

5.1. Couette flow

The analysis here is restricted to two-dimensional perturbations, $\tilde {v}^{\ast }_z = \tilde {u}^{\ast }_z = k^{\ast }_z = 0$, in order to simplify the calculation procedure. In the base state, the strain rate in the fluid, scaled by $(G/\mu _f)$, is $\varGamma = (\mu _f V_f/G h_f)$, where $V_f$ is the velocity of the top plate. The strain rate in the solid in the base state $(\textrm {d} \bar {u}^{\ast }_x/\textrm {d} Y^{\ast })$, (3.37), is also $\varGamma$. Equations (2.8)–(2.9) can be reduced to one equation for the fluid velocity by taking $(\imath k^{\ast }_x) \times$ (2.8)–$(\textrm {d} / \textrm {d} y^{\ast })$(2.9), and using the mass conservation equation (2.7) to express $\tilde {v}^{\ast }_x$ in terms of $(\textrm {d} \tilde {v}^{\ast }_y/\textrm {d} y^{\ast })$:

(5.1)\begin{equation} \left( \frac{\textrm{d}^{2}}{\textrm{d} y^{{\ast} 2}} - k^{{\ast} 2}_x \right)^{2} \tilde{v}^{{\ast}}_y = 0. \end{equation}

Equations (3.39)–(3.41) for the solid displacement field can be reduced to one equation by taking $(\imath k_x) \times$ (3.40)–$((\textrm {d}/{\textrm {d} y}) - \imath k \varGamma )$ (3.41), and using the mass conservation equation (3.39) to express $\tilde {u}^{\ast }_x$ in terms of $\tilde {u}^{\ast }_y$:

(5.2)\begin{align} & \left( (1 + \mu_r s^{{\ast}}) \left(\frac{\textrm{d}^{2}}{\textrm{d} Y^{{\ast} 2}} - k^{{\ast} 2}_x \right) - \mu_r s^{{\ast}} \left( \varGamma^{2} k^{{\ast} 2}_x + 2 \imath k^{{\ast}}_x \varGamma \frac{\textrm{d}}{\textrm{d} Y^{{\ast}}} \right) \right) \nonumber\\ &\quad \times \left( \frac{\textrm{d}^{2}}{\textrm{d} Y^{{\ast} 2}} - 2 \imath k^{{\ast}}_x \frac{\textrm{d}}{\textrm{d} Y^{{\ast}}} -k^{{\ast} 2}_x (1 + \varGamma^{2}) \right) \tilde{u}^{{\ast}}_y = 0. \end{align}

The boundary conditions, (3.43a,b)–(3.45), expressed in terms of scaled variables are

(5.3a,b)$$\begin{gather} s^{{\ast}} \tilde{u}^{{\ast}}_x = \tilde{v}^{{\ast}}_x + \varGamma \tilde{u}^{{\ast}}_y, \quad s^{{\ast}} \tilde{u}^{{\ast}}_y = \tilde{v}^{{\ast}}_y, \end{gather}$$

(5.4)$$\begin{gather} - \tilde{p}_s^{{\ast}} + 2 (1 + \mu_r s^{{\ast}} ) \frac{\textrm{d} \tilde{u}^{{\ast}}_y}{\textrm{d} Y^{{\ast}}} - 2 \mu_r \imath k^{{\ast}}_x s^{{\ast}} \tilde{u}^{{\ast}}_y \varGamma - \varGamma^{2} \tilde{u}^{{\ast}}_y ={-} \tilde{p}_f^{{\ast}} + 2 \frac{\textrm{d} \tilde{v}^{{\ast}}_y}{\textrm{d} y^{{\ast}}}, \end{gather}$$

(5.5)$$\begin{gather}( 1 + \mu_r s^{{\ast}} ) \left( \frac{\textrm{d} \tilde{u}^{{\ast}}_x}{\textrm{d} Y^{{\ast}}} + \imath k^{{\ast}}_x \tilde{u}^{{\ast}}_y \right) + \varGamma \frac{\textrm{d} \tilde{u}^{{\ast}}_y}{\textrm{d} Y^{{\ast}}} - \mu_r s^{{\ast}} \imath k^{{\ast}}_x \tilde{u}^{{\ast}}_x \varGamma = \left( \frac{\textrm{d} \tilde{v}^{{\ast}}_x}{\textrm{d} y^{{\ast}}} + \frac{\textrm{d} \tilde{v}^{{\ast}}_y}{d x^{{\ast}}} \right). \end{gather}$$

A particularly simple case is $\mu _r = 0$, where (5.2) is also independent of the growth rate $s$. Thus, both the solid and fluid equations are quasi-static, and are independent of time. The time dependence enters only through the velocity boundary conditions, (5.3a,b). In this case, the determinant of the linear stability matrix reduces to a quadratic equation for the growth rate $s$, and this can be explicitly solved to determine the growth rate. Kumaran et al. (1994) first showed that there is an instability when $\varGamma$ increases beyond a critical value even when inertia is neglected. The lowest value of $\varGamma$ for the transition from stable to unstable modes is at a finite value of the wavenumber $k^{\ast }_x$, as shown in figure 6(a) for the case where $\mu _r = 0$, and so this instability is a finite-wavenumber instability.

Figure 6. (a) The real part of the growth rate of the most unstable mode as a function of the wavenumber $k^{\ast }_x$ for the viscous instability for $H=1$, $\mu _r = 0$ and $\varGamma = 0$ ($\circ$), $0.5$ ($\triangle$), $0.8$ ($\boldsymbol {\nabla }$), $0.9$ ($\Box$) and $1.0$ ($\diamond$). (b) Neutral stability curves plotted in the $(\mu _r/H^{2})$ versus $(\varGamma _c H)$ parameter space for $H=1$ ($\circ$), $10$ ($\triangle$), $50$ ($\boldsymbol {\nabla }$), $100$ ($\Box$) and 500 $(\diamond )$ for the viscous instability in the flow past a neo-Hookean solid. The dashed curve in (b) is the result for the linear model. Panel (b) is redrawn, with permission, from Chokshi & Kumaran (2008a).

There is an instability even for $\mu _r > 0$, and this instability is caused due to the second term on the right-hand side of (5.3a,b) for the tangential velocity, which is the change in the mean velocity of the interface due to interface displacement. There is a discontinuity in the strain rate at the interface, in a manner similar to the discontinuity in the velocity at the interface which causes the Kelvin–Helmholtz instability. The strain rate in the fluid is non-zero at the interface, while the strain rate in the solid is zero because the solid is stationary. When the interface is displaced, the mean velocity on the fluid side at the interface is non-zero, whereas the mean velocity on the solid side is zero. This results in the second term on the right-hand side in the tangential velocity boundary condition, (5.3a,b), which destabilises the flow when the parameter $\varGamma$ exceeds a critical value. The instability of mechanism is the transfer of energy from the mean flow to the fluctuations due to the shear work done at the interface. Figure 6(a) shows the results for the real part of the scaled growth rate $s^{\ast }$ as a function of the wavenumber $k^{\ast }_x$ for different values of the parameter $\varGamma$. The flow is always stable in the long-wave limit, $k^{\ast } \ll 1$, but there is an instability at a finite wavenumber when $\varGamma$ exceeds a critical value.

The analysis in Kumaran et al. (1994) employed the linear model for the solid constitutive relation, with the nonlinear term in the definition of the strain tensor (last term on the right-hand side in (3.11)) neglected. The difficulty with the linear model is that the constitutive relation is not invariant under rotation. Gkanis & Kumar (2003) first used the neo-Hookean model, which is rotationally invariant strain ((3.11) including the nonlinear term), in their stability analysis in the viscous limit, and included the effect of the normal stress difference in the base state (the difference in the diagonal terms in (3.32), which results in the fourth term on the left-hand side in the tangential stress balance (5.4)). The qualitative nature of the instability due to the shear work at the interface is the same for the neo-Hookean model, though the critical strain rate for the instability is lower than that for the linear model. The difference in the results for the linear and neo-Hookean models decreased to zero in the limit $H \gg 1$, where the thickness of the wall is much smaller than the fluid thickness. The normal stress difference caused an additional short-wave instability in the limit of zero Reynolds number, which is damped by the effect of surface tension. Chokshi & Kumaran (2008a) also analysed the stability of the flow of a viscous fluid past a neo-Hookean solid at zero Reynolds number with a rotationally invariant strain, and found that the critical strain rate for the neo-Hookean model is lower than that for the linear model. The linear model (Kumaran et al. 1994) predicted that there is an instability only for $H > \sqrt {\mu _r}$, and the flow is always stable for $H < \sqrt {\mu _r}$. For the neo-Hookean solid model, Chokshi & Kumaran (2008a) found that there is an instability for all values of $\mu _r$, and an increase in the solid viscosity $\mu _r$ does not monotonically increase the critical value $\varGamma _c$. An example of the result from Chokshi & Kumaran (2008a) is shown in figure 6(b), where $\varGamma _c H$ is plotted as a function of $(\mu _r/H^{2})$. Whereas the value of $\varGamma _c$ diverges for the linear model at $(\mu _r/H^{2}) = 1$, it is always finite for the neo-Hookean model.

Chokshi & Kumaran (2008a) also carried out a weakly nonlinear analysis to determine whether the bifurcation is supercritical or subcritical, using the neo-Hookean model for the solid. Both the Eulerian and Lagrangian formulations were used, and it was shown that the results of the two approaches are in quantitative agreement. An earlier weakly nonlinear analysis by Shankar & Kumaran (2001b), which employed the linear strain measure ((3.11) without the last term on the right-hand side), arrived at qualitatively similar results. The regions in the $H$–$\mu _r$ parameter space for the supercritical and subcritical instabilities were mapped out in Chokshi & Kumaran (2008a), and it was found that the instability is subcritical in most of the parameter space, though there are isolated regions where the instability is supercritical. However, the reduction in $\varGamma _c$ upon inclusion of the second- and third-order terms in the amplitude equation is relatively small, indicating that the linear analysis provides an accurate estimate of $\varGamma _c$. The theoretical predictions for the value of $\varGamma _c$ and the nature of the bifurcation are in agreement with experimental results; the comparison is discussed in § 8. Thus, the viscous instability in the Couette flow past a compliant surface is now a well-understood phenomenon.

Though the stability of a viscous flow past a viscoelastic wall has been relatively well studied, there are fewer studies of other wall models. Thaokar et al. (2001) used (3.1) and (3.2) discussed in § 3.1, and showed that tangential motion of the interface has a qualitative effect on the flow stability at zero Reynolds number. The stability of the flow past an infinitesimal membrane in the low-Reynolds-number limit was considered by Thaokar & Kumaran (2002), using the model discussed in § 3.2. This configuration was found to be unstable for all non-zero values of the fluid velocity due to a long-wavelength instability – the growth rate of perturbations increased proportional to $k^{2}$ when the wavenumber $k$ is small.

5.2. Pressure-driven flow

The stability of a pressure-driven flow is found to be extremely sensitive to the model used for the solid and the nature of the boundary conditions. In the absence of inertia, the fluid equation is the same as that for a Couette flow, (5.1), but the solid equation cannot be reduced to a fourth-order equation with constant coefficients.

The first study of the pressure-driven flow in a pipe (Kumaran 1995) employed the linear strain measure ((3.11) without the last term on the right-hand side) for the solid constitutive relation which is not invariant under rotation. A linear instability was predicted when the parameter $\varGamma$ exceeded a critical value. The mechanism, the transport of energy from the mean flow to the fluctuations due to the shear work at the interface, is the same as that for a Couette flow. The neo-Hookean wall model was used by Gkanis & Kumar (2005), with strain given by (3.11), for the pressure-driven flow in a channel lined with one or two compliant walls, as well as a combined Couette/Poiseuille flow. They reported that the Couette flow is easier to destabilise than the Poiseuille flow, but the Poiseuille flow does become unstable when the scaled characteristic velocity exceeds a critical value. The normal stress difference in the solid in the base state causes a short-wave instability, which is stabilised by surface tension. At finite wavenumber, there is an instability due to the shear work done at the interface. In addition, there is a term in the normal stress boundary condition due to the variation in the normal stress with height in the base state. This term arises because the boundary conditions were applied at the perturbed surface on the solid side, and not on the unperturbed surface as appropriate for the Lagrangian formulation. This additional term in the normal stress equation has a stabilising effect, and due to this, the Poiseuille flow is more stable than the Couette flow.

The flow in a deformable channel in the viscous limit was also studied by Gaurav & Shankar (2010) using the rotational invariant strain measure (equation (3.11)), who reported that the viscous flow in a channel is unstable only to the short-wavelength modes caused by the normal stress difference. They did not find an instability at finite wavenumber due to the shear work done at the interface. Gaurav (2009) came to the same conclusion for the viscous flow in a tube which has an annular wall made of a neo-Hookean solid. The models used by Gkanis & Kumar (2005) and Gaurav & Shankar (2010) are exactly the same, but there is a difference in the continuity conditions at the solid–fluid interface. On the fluid side, the variation of the shear stress with cross-stream distance, which is the second term on the right-hand side in (3.45), was not included in Gkanis & Kumar (2005), but was included in Gaurav & Shankar (2010). On the solid side, Gkanis & Kumar (2005) included a term due to the gradient of the normal stress difference in the base state in the normal stress balance equation; this was not included in Gaurav & Shankar (2010) because the boundary conditions were applied at the undeformed interface in the Lagrangian formulation. Neither Gkanis & Kumar (2005) nor Gaurav & Shankar (2010) included the streamwise contribution to the normal stress balance due to the pressure gradient, which is the last term on the right-hand side in (3.44). The linear stability analysis of Gaurav & Shankar (2010) and Gaurav (2009) was based on a Taylor expansion of the quantities about the equilibrium state (figure 4a) and not about the strained base state (figure 4b). These issues were revisited in Patne et al. (2017) and Patne & Shankar (2019a), who used a consistent formulation where the perturbation is applied on the base steady state. The qualitative nature of the instability does not change, and the high-wavenumber modes due to the normal stress difference are the critical modes in this case.

The studies conducted so far reveal that the stability of a pressure-driven flow in a channel or tube with compliant walls is extremely sensitive to the wall model used, and small differences in the interface conditions qualitatively change the stability characteristics. The stability of the pressure-driven flow is also difficult to verify experimentally, since the parameter regime of low Reynolds number and finite $\varGamma$ requires a very large pressure difference for driving the flow, and soft materials will not be able to sustain the pressure generated. This is discussed in § 8.1.

6. High-Reynolds-number instability

There are three kinds of instabilities for a flow past a compliant wall at high Reynolds number. The first are the ‘inviscid’ modes, where the instability is predicted by the inviscid Euler equations when the viscous effects are neglected. The second are structurally similar to the Tollmien–Schlichting modes which destabilise the flow in a rigid channel. These involve an internal critical layer of thickness ${Re}^{-1/3}$ within the flow at a location where the flow velocity is equal to the wave velocity. At high Reynolds number, the flow becomes unstable due to the viscous effects in the critical layer. The third are the wall modes. These are present for the flow in a rigid channel/tube (Gill 1965), but are always stable. For the wall modes, the effect of viscosity is significant in a layer of thickness ${Re}^{-1/3}$ at the wall of the channel. This is in contrast to the Tollmien–Schlichting modes, where the effect of viscosity is important in an internal critical layer of thickness ${Re}^{-1/3}$ within the flow. The mechanism of destabilisation of the wall modes is the transport of energy from the mean flow to the fluctuations due to the shear work done at the interface; this mechanism is present only when the wall is deformable. In this respect, the wall mode instability is qualitatively different from, and is not a simple continuation of, the Tollmien–Schlichting instability in a rigid tube/channel.

6.1. Theorems for inviscid flow

Equivalents of the classical theorems for inviscid flows (Drazin & Reid 1981) can be derived for the spring-backed wall model, but with restrictions, as indicated in table 2. Inertial scaling is used for the fluid velocity and pressure in the present section. The scaled fluid and wave velocity are defined as $\tilde {v}_i^{{{\dagger}} } = (\tilde {v}_i/V_f)$, $c^{{{\dagger}} } = (c/V_f)$ and $t^{{{\dagger}} } = (t V_f / h_f)$, where $V_f$ is the characteristic velocity. The fluid pressure is defined as $p_f^{{{\dagger}} } = p_f / (\rho _f V_f^{2})$. All other quantities are non-dimensionalised as shown in table 1.

Table 2. The conditions for the existence of unstable modes in channel and tube flows when the mean velocity is zero at the compliant walls, and all walls with non-zero mean velocity are rigid. The function ${\mathcal {G}} =\frac {r^{\ast }}{n^{2} + k^{\ast 2}_x r^{\ast 2}} \frac {\textrm {d} \bar {v}_x^{{{\dagger}} }}{\textrm {d} r^{\ast }}$, the velocity $\bar {v}_{xL}^{{{\dagger}} } = \mbox {Min}(\mbox {Min}(\bar {v}_x^{{{\dagger}} }), 0)$ and $\bar {v}_{xU}^{{{\dagger}} } = \mbox {Max}(\mbox {Max}(\bar {v}_x^{{{\dagger}} }), 0)$.

The wall displacement $u^{\ast }$ is zero in the base state. The wall displacement due to the fluid perturbation is expressed in a manner similar to that for the fluid in (2.6):

(6.1)\begin{equation} u^{{\ast}} = \tilde{u}^{{\ast}} \exp{(\imath k^{{\ast}}_x (x^{{\ast}} - c^{{{\dagger}}} t^{{{\dagger}}}) + \imath k^{{\ast}}_z z^{{\ast}})}. \end{equation}

The boundary conditions at the surface are the continuity of normal velocity and the continuity of normal stress at the perturbed surface. When these conditions are expressed using a Taylor expansion about their values at the flat wall, the normal velocity condition is

(6.2)\begin{equation} \tilde{v}_y^{{{\dagger}}} ={-} \imath k^{{\ast}}_x c^{{{\dagger}}} \tilde{u}^{{\ast}}. \end{equation}

The relation (3.1) between the fluid pressure and the wall displacement is combined with the above normal velocity continuity (6.2) to provide a relation between the pressure and the normal velocity at the surface:

(6.3)\begin{equation} \tilde{p}_f^{{{\dagger}}} ={\pm} \frac{\imath (E^{{\ast}} - k^{{\ast} 2}_x c^{{{\dagger}} 2} I^{{\ast}} - \imath k^{{\ast}}_x c_x^{{{\dagger}}} D^{{\ast}} + T^{{\ast}} (k^{{\ast} 2}_x + k^{{\ast} 2}_z)) \tilde{v}_y^{{{\dagger}}}}{k^{{\ast}}_x c^{{{\dagger}}}}, \end{equation}

where the positive sign is applicable at the upper boundary $y^{\ast } = y^{\ast }_u$ in figure 3(b) where the unit normal is in the $-y$ direction and the negative sign is applicable at the lower boundary $y^{\ast } = y^{\ast }_l$ where the unit normal is in the $+y$ direction.

A brief derivation of some of the classical theorems is given here for the flow in a channel with compliant walls; these and the equivalent results for a tube flow are summarised in table 2. The objective is to demonstrate how the wall dynamics is incorporated in the calculation, and to show how the restrictions arise on the wall dynamics. The important restriction for the channel flow is that either the mean velocity should be zero at the surface if it is compliant or the surface should be rigid if the mean velocity is non-zero. Thus, the present analysis will apply to configurations of the type shown in figure 1, but the theorems would not apply in general when there are two compliant surfaces moving relative to each other. In the case of a tube flow, the important restriction is that the mean velocity at the tube surface should be zero. The detailed calculation for a rigid channel is given in Drazin & Reid (1981) and the adaptation for a channel with a compliant wall is given in Yeo & Dowling (1987). The derivation of the classical theorems for a rigid tube was first carried out by Howard & Gupta (1962) and Maslowe (1974), and the calculation for a tube with a compliant wall is given in Kumaran (1996) and Shankar & Kumaran (2000).

To derive the equivalent of Squire's theorem, add the $x$ momentum equation (2.8) and $(k^{\ast }_z / k^{\ast }_x) \times$ the $z$ momentum equation (2.10), and use the transformation

(6.4a–d)\begin{equation} k_x^{\unicode{x2021}} = \sqrt{k^{{\ast} 2}_x + k^{{\ast} 2}_z}, \quad k_x^{\unicode{x2021}} \tilde{v}_x^{\unicode{x2021}} = k^{{\ast}}_x \tilde{v}_x^{{{\dagger}}} + k^{{\ast}}_z \tilde{v}_z^{{{\dagger}}}, \quad \tilde{p}_f^{\unicode{x2021}} = (\tilde{p}^{{{\dagger}}}_f k^{{\ast}}_x / k_x^{\unicode{x2021}}), \quad {Re}^{\unicode{x2021}} = {Re} (k^{{\ast}}_x/k_x^{\unicode{x2021}})\end{equation}

to obtain the transformed $x$ momentum conservation equation:

(6.5)\begin{equation} \imath k_x^{\unicode{x2021}} (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}}) \tilde{v}_x^{\unicode{x2021}} + \tilde{v}_y^{{{\dagger}}} \frac{\textrm{d} \bar{v}^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} ={-} \imath k_x^{\unicode{x2021}} \tilde{p}_f^{\unicode{x2021}} + \frac{1}{{Re}^{\unicode{x2021}}} \left(\frac{\textrm{d}}{\textrm{d} y^{{\ast} 2}} - k^{\unicode{x2021} 2} \right) \tilde{v}_x^{\unicode{x2021}}.\end{equation}

The transformation (6.4a–d) is applied on the mass conservation equation (2.7) to obtain

(6.6)\begin{equation} \imath k_x^{\unicode{x2021}} \tilde{v}_x^{\unicode{x2021}} + \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} = 0. \end{equation}

The $y$ momentum conservation equation, (2.9), is multiplied by $(k_x^{\unicode{x2021} }/ k^{\ast }_x)$ to obtain

(6.7)\begin{equation} \imath k_x^{\unicode{x2021}} (\bar{v}_x^{{{\dagger}}}(y) - c^{{{\dagger}}}) \tilde{v}_y^{{{\dagger}}} ={-} \frac{\textrm{d} \tilde{p}_f^{\unicode{x2021}}}{\textrm{d} y^{{\ast}}} + \frac{1}{{Re}^{\unicode{x2021}}}\left( \frac{\textrm{d}^{2}}{\textrm{d} y^{{\ast} 2}} - k^{\unicode{x2021} 2}_x \right) \tilde{v}_y^{{{\dagger}}}.\end{equation}

Equations (6.5)–(6.7) are the linear stability equations for a two-dimensional flow with wavenumber $k_x^{\unicode{x2021} }$; from these, it can be inferred that if a three-dimensional perturbation with wavenumber $(k^{\ast }_x,k^{\ast }_z)$ has growth rate $k^{\ast }_x c^{{{\dagger}} }$ when the Reynolds number is ${Re}$, a two-dimensional perturbation with wavenumber $k_x^{\unicode{x2021} }$ has a higher growth rate $k_x^{\unicode{x2021} } c^{{{\dagger}} }$ at a lower Reynolds number ${Re}^{\unicode{x2021} } = ({Re} k^{\ast }_x/ k_x^{\unicode{x2021} })$. This leads to Squire's theorem, which states that two-dimensional disturbances are always more unstable than three-dimensional ones, and therefore the stability limit can be identified by examining two-dimensional disturbances alone. For a flow bounded by a spring-backed plate, it is necessary to transform the parameters in (6.3) as well,

(6.8a,b)\begin{equation} E^{\unicode{x2021}} = (E^{{\ast}} k^{\unicode{x2021} 2}_x / k^{{\ast} 2}_x), \quad D^{\unicode{x2021}} = (D^{{\ast}} k_x^{\unicode{x2021}}/k^{{\ast}}_x), \end{equation}

to obtain the condition equivalent to (6.3) for the two-dimensional perturbation with wavenumber $k_x^{\unicode{x2021} }$:

(6.9)\begin{equation} \tilde{p}^{\unicode{x2021}} ={\pm} \frac{\imath (E^{\unicode{x2021}} - k^{\unicode{x2021} 2}_x c^{{{\dagger}} 2} I^{{\ast}} - \imath k_x^{\unicode{x2021}} c_x^{{{\dagger}}} D^{{{\dagger}}} + T^{\ast}k^{\unicode{x2021} 2}) \tilde{v}_y^{{{\dagger}}}}{k_x^{\unicode{x2021}} c^{{{\dagger}}}}. \end{equation}

Therefore, for a parallel shear flow with Reynolds number ${Re}$ and wall parameters $E^{\ast }, I^{\ast }, D^{\ast }$ and $T^{\ast }$, a three-dimensional disturbance with wavenumbers $k^{\ast }_x, k^{\ast }_z$ has growth rate $- \imath k^{\ast }_x c^{{{\dagger}} }$; a two-dimensional disturbance with wavenumber $k_x^{\unicode{x2021} }$ and lower Reynolds number ${Re}^{\unicode{x2021} } = ({Re} k^{\ast }_x/ k_x^{\unicode{x2021} })$ with higher wall elasticity and dissipation, $(E^{\ast } k^{\unicode{x2021} 2}_x / k^{\ast 2}_x), (D^{\ast } k_x^{\unicode{x2021} } / k^{\ast }_x)$ and the same inertia parameter $I^{\ast }$ and surface tension $T^{\ast }$, has a growth rate $- \imath k_x^{\unicode{x2021} } c^{{{\dagger}} }$ that is larger in magnitude. Thus, Squire's theorem relates the stability of a three-dimensional flow to a two-dimensional flow with lower Reynolds number and higher wall elasticity and dissipation.

There is no equivalent of Squire's theorem for the flow in a pipe with either rigid or compliant walls.

For two-dimensional disturbances, the Orr–Sommerfeld equation is obtained by adding $(\textrm {d} /\textrm {d} y^{\ast }) \times$ (2.8) and $- \imath k^{\ast }_x \times$ (2.9), and using the mass conservation equation (2.7) to express $\tilde {v}_x^{{{\dagger}} }$ in terms of $\tilde {v}_y^{{{\dagger}} }$:

(6.10)\begin{equation} \imath k^{{\ast}}_x (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}}) \left( \frac{\textrm{d}^{2}}{\textrm{d} y^{{\ast} 2}} - k^{{\ast} 2}_x \right) \tilde{v}_x^{{{\dagger}}} - \imath k^{{\ast}}_x \tilde{v}_y^{{{\dagger}}} \frac{\textrm{d}^{2} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast} 2}} = \frac{1}{{Re}} \left( \frac{\textrm{d}^{2}}{\textrm{d} y^{{\ast} 2}} - k^{{\ast} 2}_x \right)^{2} \tilde{v}_y^{{{\dagger}}}.\end{equation}

The Orr–Sommerfeld equation is usually expressed in terms of the stream function, but here it is expressed in terms of the cross-stream velocity fluctuation so that it is easier to apply the boundary conditions. An inviscid flow, where viscous effects are neglected, is described by the Rayleigh equation:

(6.11)\begin{equation} \left( \frac{\textrm{d}^{2}}{\textrm{d} y^{{\ast} 2}} - k^{{\ast} 2}_x \right) \tilde{v}_y^{{{\dagger}}} - \frac{\tilde{v}_y^{{{\dagger}}}}{\bar{v}^{{{\dagger}}} - c^{{{\dagger}}}} \frac{\textrm{d}^{2} \bar{v}^{{{\dagger}}}}{\textrm{d} y^{{\ast} 2}} = 0.\end{equation}

The pressure is related to the velocity through $x$ momentum equation (2.8) in which the viscous term is neglected:

(6.12)\begin{equation} \tilde{p}^{{{\dagger}}} = \frac{(\bar{v}^{{{\dagger}}} - c^{{{\dagger}}})}{\imath k_x^{{\ast}}} \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} - \frac{\tilde{v}_y^{{{\dagger}}}}{\imath k^{{\ast}}_x} \frac{\textrm{d} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast}}}.\end{equation}

At the compliant boundary, the pressure is related to the normal velocity $\tilde {v}^{{{\dagger}} }$ by (6.3). Substituting for the pressure from (2.8) into (6.12), we obtain a relation between the normal velocity and its derivative at the compliant boundary:

(6.13)\begin{equation} \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} = \frac{\tilde{v}_y^{{{\dagger}}}}{\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}}} \frac{\textrm{d} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} \mp \frac{(E^{{\ast}} - k^{{\ast} 2}_x c^{{{\dagger}} 2} I^{{\ast}} - \imath k^{{\ast}}_x c^{{{\dagger}}} D^{{\ast}} + T^{{\ast}} k^{{\ast} 2}_x) \tilde{v}_y^{{{\dagger}}}}{c^{{{\dagger}}} (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}})}, \end{equation}

where the negative (positive) sign multiplies the second term on the right-hand side at the upper (lower) boundary in figure 3(b).

The theorems derived here are applicable when the mean velocity at the wall is zero in the base state when the surface is compliant, as shown in figure 1(b), or when the moving boundary is rigid, as shown in figure 1(a) for the top moving surface. In the former case, the mean velocity in (6.13) is equal to zero at the boundaries, while in the latter $\tilde {v}_x^{{{\dagger}} }$ and $\tilde {v}_y^{{{\dagger}} }$ are both zero at the boundary. Therefore, (6.13) is simplified by setting the mean velocity equal to zero in the following analysis.

If we multiply the Rayleigh equation (6.11) by $\tilde {v}^{{{\dagger}} cc}_y$, the complex conjugate of $\tilde {v}_y^{{{\dagger}} }$, and integrate across the channel, we obtain

(6.14)\begin{equation} \left. \tilde{v}^{{{\dagger}} cc}_y \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y} \right|_{y^{{\ast}}=y^{{\ast}}_l}^{y^{{\ast}}=y^{{\ast}}_u} - \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} \left( \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} \frac{\textrm{d} \tilde{v}^{{{\dagger}} cc}_y}{\textrm{d} y^{{\ast}}} + \frac{\tilde{v}_y^{{{\dagger}}} \tilde{v}^{{{\dagger}} cc}_y}{\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}}} \frac{\textrm{d}^{2} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast} 2}} + k^{{\ast} 2}_x \tilde{v}_y^{{{\dagger}}} \tilde{v}^{{{\dagger}} cc}_y \right) = 0. \end{equation}

Substituting the value of $(\textrm {d} \tilde {v}^{{{\dagger}} }/\textrm {d} y^{\ast })$ at the boundary from (6.13), and setting $\bar {v}_x = 0$ at the boundaries, we obtain

(6.15)\begin{align} &- \left. \frac{\tilde{v}_y^{{{\dagger}}} \tilde{v}^{{{\dagger}} cc}_y}{c^{{{\dagger}}}} \frac{\textrm{d} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} \right|_{y^{{\ast}}_l}^{y^{{\ast}}_u} + \sum_{y^{{\ast}}_l,y^{{\ast}}_u} \frac{(E^{{\ast}} + T^{{\ast}} k^{{\ast} 2}_x - k^{{\ast} 2}_x c^{{{\dagger}} 2} I^{{\ast}} - \imath k^{{\ast}}_x c^{{{\dagger}}} D^{{\ast}}) \tilde{v}_y^{{{\dagger}}} \tilde{v}^{{{\dagger}} cc}_y}{c^{{{\dagger}} 2}} \nonumber\\ &\quad - \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} \left( \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} \frac{\textrm{d} \tilde{v}^{{{\dagger}} cc}_y}{\textrm{d} y^{{\ast}}} + \frac{\tilde{v}_y^{{{\dagger}}} \tilde{v}^{{{\dagger}} cc}_y}{\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}}} \frac{\textrm{d}^{2} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast} 2}} + k^{{\ast} 2}_x \tilde{v}_y^{{{\dagger}}} \tilde{v}^{{{\dagger}} cc}_y\right) = 0. \end{align}

The Rayleigh theorem is derived by multiplying (6.15) by $c^{{{\dagger}} }$, and taking the imaginary part of the resulting equation, to obtain

(6.16)\begin{align} &- c^{{{\dagger}}}_I \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} \left( \left| \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} \right|^{2} + \frac{|\tilde{v}_y^{{{\dagger}}}|^{2} \bar{v}_x^{{{\dagger}}}}{|\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}}|^{2}} \frac{\textrm{d}^{2} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast} 2}} + k^{{\ast} 2}_x |\tilde{v}_y^{{{\dagger}}}|^{2} \right) \nonumber\\ &\quad - \sum_{y^{{\ast}}_l,y^{{\ast}}_u} \left[ c^{{{\dagger}}}_I \left(\frac{(E^{{\ast}} + T^{{\ast}} k^{{\ast} 2}_x)}{|c^{{{\dagger}}}|^{2}} + I^{{\ast}} k^{{\ast} 2}_x \right) + k^{{\ast}}_x D^{{\ast}} \right] |\tilde{v}_y^{{{\dagger}}}|^{2} = 0. \end{align}

The second term on the left-hand side in the above equation is negative for unstable modes with $c^{{{\dagger}} }_I > 0$, because the term in the square bracket is positive. Therefore, unstable modes with $c^{{{\dagger}} }_I > 0$ can exist only if the first term on the left-hand side is positive, that is, if $\bar {v}_x^{{{\dagger}} } (\textrm {d}^{2} \bar {v}_x^{{{\dagger}} }/\textrm {d} y^{\ast 2})$ is negative somewhere in the flow. This is the equivalent of the Rayleigh theorem for the flow past a spring-backed plate. A corollary is that if $\bar {v}_x^{{{\dagger}} } (\textrm {d}^{2} \bar {v}_x^{{{\dagger}} }/\textrm {d} y^{\ast 2}) > 0$ everywhere in the domain, then perturbations are neutrally stable ($c^{{{\dagger}} }_I = 0$) if $D^{\ast }$ is zero, and stable $(c^{{{\dagger}} }_I < 0)$ if $D^{\ast }$ is positive.

The Fjørtoft theorem is derived from the imaginary part of (6.15):

(6.17)\begin{align} &\left. \frac{|\tilde{v}_y^{{{\dagger}}}|^{2} c^{{{\dagger}}}_I}{|c^{{{\dagger}}}|^{2}} \frac{\textrm{d} \bar{v}_y^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} \right|_{y^{{\ast}}_l}^{y^{{\ast}}_u} - c^{{{\dagger}}}_I \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} \frac{|\tilde{v}_y^{{{\dagger}}}|^{2}}{|\bar{v}^{{{\dagger}}} - c^{{{\dagger}}}|^{2}} \frac{\textrm{d}^{2} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast} 2}} \nonumber\\ &\quad - \sum_{y^{{\ast}}_l,y^{{\ast}}_u} \frac{c_R^{{{\dagger}}}}{|c^{{{\dagger}}}|^{2}} \left( \frac{2 c^{{{\dagger}}}_I (E^{{\ast}} + T^{{\ast}} k^{{\ast} 2}_x)}{|c^{{{\dagger}}}|^{2}} + k^{{\ast}}_x D^{{\ast}} \right) |\tilde{v}_y^{{{\dagger}}}|^{2} = 0. \end{align}

For a downstream/upstream travelling wave, $c_R^{{{\dagger}} } \gtrless 0$, if the first term on the left-hand side is negative/positive, an unstable mode with $c^{{{\dagger}} }_I > 0$ can exist only if $(\textrm {d}^{2} \bar {v}_x^{{{\dagger}} }/\textrm {d} y^{\ast 2}) \lessgtr 0$ somewhere in the flow. This leads to a rather restrictive statement for the Fjørtoft theorem in table 2.

A bound on the real part of $c_R$ can be derived by defining the function $\tilde {g} = (\tilde {v}_y^{{{\dagger}} } / (\bar {v}_x^{{{\dagger}} } - c))$. The Rayleigh equation (6.11) is be expressed in terms of $\tilde {g}$, multiplied by $(\bar {v}_x^{{{\dagger}} } - c^{{{\dagger}} })$ and simplified:

(6.18)\begin{equation} \frac{\textrm{d}}{\textrm{d} y^{{\ast}}} \left( (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}})^{2} \frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}} \right) - k^{{\ast} 2}_x (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}})^{2} \tilde{g} = 0. \end{equation}

The boundary condition (6.13) can be expressed in terms of $\tilde {g}$:

(6.19)\begin{equation} \frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}} ={\pm} \frac{(E^{{\ast}} + T^{{\ast}} k^{{\ast} 2}_x - k^{{\ast} 2}_x c^{{{\dagger}} 2} I^{{\ast}} - \imath k^{{\ast}}_x c^{{{\dagger}}} D^{{\ast}}) \tilde{g}}{c^{{{\dagger}} 2}}. \end{equation}

Equation (6.19) is is multiplied by $\tilde {g}^{cc}$, the complex conjugate of $\tilde {g}$, and integrated over the height of the channel to obtain

(6.20)\begin{equation} \left. \tilde{g}^{cc} c^{{{\dagger}} 2} \frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}} \right|_{y^{{\ast}}=y^{{\ast}}_l}^{y^{{\ast}} = y^{{\ast}}_u} - \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}})^{2} \left( \frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}} \frac{\textrm{d} \tilde{g}^{cc}}{\textrm{d} y^{{\ast}}} + k^{{\ast} 2}_x \tilde{g} \tilde{g}^{cc} \right) = 0. \end{equation}

Here, the mean velocity is assumed to be zero at the boundaries. Substituting (6.19) for $(\textrm {d} \tilde {g} / \textrm {d} y^{\ast })$ at the boundaries, we obtain

(6.21)\begin{equation} \sum_{y^{{\ast}}_l,y^{{\ast}}_u} (E^{{\ast}} - k^{{\ast} 2}_x c^{{{\dagger}} 2} I^{{\ast}} - \imath k^{{\ast}}_x c^{{{\dagger}}} D^{{\ast}}) |\tilde{g}|^{2} - \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} (\bar{v}_x^{{{\dagger}}} - c^{{{\dagger}}})^{2} \left( \left| \frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}} \right|^{2} + k^{{\ast} 2}_x |\tilde{g}|^{2} \right) = 0.\end{equation}

The imaginary part of (6.21) is

(6.22)\begin{equation} 2 c^{{{\dagger}}}_I \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} (\bar{v}_x^{{{\dagger}}} - c_R^{{{\dagger}}} ) \left( \left| \frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}} \right|^{2} + k^{{\ast} 2}_x |\tilde{g}|^{2} \right) = \sum_{y^{{\ast}}_l,y^{{\ast}}_u} c_R^{{{\dagger}}} (2 k^{{\ast} 2}_x c^{{{\dagger}}}_I I^{{\ast}} + k^{{\ast}}_x D^{{\ast}} ) |\tilde{g}|^{2}.\end{equation}

For unstable modes with $c^{{{\dagger}} }_I > 0$, the right-hand side of the above equation has the same sign as $c_R^{{{\dagger}} }$. For $c_R^{{{\dagger}} } > 0$, the right-hand side is positive, and a necessary condition for the existence of solutions is $\mbox {Max}(\bar {v}_x^{{{\dagger}} }) > 0$ and $c_R^{{{\dagger}} } < \mbox {Max}(\bar {v}_x^{{{\dagger}} })$. For $c_R^{{{\dagger}} } < 0$, the right-hand side is negative, and a necessary condition for the existence of solutions is that $\mbox {Min}(\bar {v}_x^{{{\dagger}} }) < 0$ and $c_R^{{{\dagger}} } > \mbox {Min}(\bar {v}_x^{{{\dagger}} })$. From this, we obtain the condition in the seventh row of table 2, that $\bar {v}_{xL}^{{{\dagger}} } < c_R^{{{\dagger}} } < \bar {v}_{xU}^{{{\dagger}} }$, where $\bar {v}_{xL}^{{{\dagger}} } = \mbox {Min}(\mbox {Min}(\bar {v}_x^{{{\dagger}} }),0)$ and $\bar {v}_{xU}^{{{\dagger}} } = \mbox {Max}(\mbox {Max}(\bar {v}_x^{{{\dagger}} }), 0)$.

Equation (6.21) is multiplied by $c^{{{\dagger}} cc}$, the complex conjugate of $c^{{{\dagger}} }$, and the imaginary part of the resulting equation is

(6.23)\begin{align} &- c^{{{\dagger}}}_I \int_{y^{{\ast}}_l}^{y^{{\ast}}_u} \, \textrm{d} y^{{\ast}} (|c^{{{\dagger}}}|^{2} - \bar{v}_{x}^{{{\dagger}} 2}) \left( \left|\frac{\textrm{d} \tilde{g}}{\textrm{d} y^{{\ast}}}\right|^{2} + k^{{\ast} 2}_x |\tilde{g}|^{2} \right) \nonumber\\ &\quad = \sum_{y^{{\ast}}_l,y^{{\ast}}_u} ( (E^{{\ast}} + T^{{\ast}} k^{{\ast} 2}_x) c^{{{\dagger}}}_I + k^{{\ast} 2}_x |c^{{{\dagger}}}|^{2} c^{{{\dagger}}}_I I^{{\ast}} + k^{{\ast}}_x |c^{{{\dagger}}}|^{2} D^{{\ast}}) |\tilde{g}|^{2} ). \end{align}

The right-hand side of the above equation is positive for unstable modes with $(c^{{{\dagger}} }_I > 0)$, and therefore a necessary condition for the existence of such modes is that $|c^{{{\dagger}} }|^{2} < \mbox {Max}(\bar {v}_{x}^{{{\dagger}} 2})$. This provides the condition in the eighth row of table 2, that $|c^{{{\dagger}} }|^{2} < \mbox {Max}(\bar {v}_{x}^{{{\dagger}} 2})$ somewhere in the flow.

For a channel with compliant walls, the Rayleigh theorem requires that $\bar {v}_x^{{{\dagger}} } (\textrm {d}^{2} \bar {v}_x^{{{\dagger}} }/\textrm {d} y^{\ast 2})$ is negative somewhere in the flow; this is in contrast to the flow in a rigid channel where $(\textrm {d}^{2} \bar {v}_x^{{{\dagger}} } / \textrm {d} y^{\ast 2})$ passes through zero somewhere in the flow. Therefore, this is not an ‘inflection point’ theorem. The parabolic flow in a channel could be unstable in the inviscid limit because the curvature is negative, but the Couette flow is always stable. The Fjørtoft theorem is much more restrictive in comparison with that for a rigid channel, since there are constraints on the slope of the velocity profiles at the boundaries. The bounds on the real part and the magnitude of the wave speed for unstable modes, and the Howard and Høiland theorems (Yeo & Dowling 1987) are identical to those for the flow in a rigid channel. These are useful, because the range of the wavenumber for unstable solutions is restricted.

The equivalent of the classical theorems for the flow in a tube with spring-backed wall are provided in table 2. The reader is referred to Shankar & Kumaran (2000) for the derivation. The practical implications of these theorems are as follows. The parabolic flow in a cylindrical tube is stable to axisymmetric disturbances, because $(\textrm {d} {\mathcal {G}} / \textrm {d} r^{\ast })$ is zero for a parabolic flow with $n=0$. However, the parabolic flow could be unstable to non-axisymmetric disturbances. The entrance flow in a tube and the flow in a converging tube could be unstable to axisymmetric and non-axisymmetric disturbances, since these have a lower curvature at the centre and higher gradients at the walls, resulting in a negative value for $(\textrm {d} {\mathcal {G}} / \textrm {d} r^{\ast })$. However, the flow in a diverging tube is always stable, because it turns out that $(\textrm {d} {\mathcal {G}} / \textrm {d} r^{\ast })$ is positive. Though these results are obtained only for a simple wall model, subsequent studies show that the results appear to hold even when the calculations are carried out for a viscoelastic continuum wall.

6.2. Internal critical layer

The limits for the magnitude, real and imaginary parts of the wave speed in the seventh to tenth rows of table 2 significantly simplify the search for solutions for the growth rates. However, these introduce a complication for neutral modes with $c^{{{\dagger}} }_I = 0$. From the condition on the seventh row of table 2, the wave velocity is in the interval $0 \leq c_R^{{{\dagger}} } \leq \mbox {Max}( \bar {v}_x)$. In this case, the Rayleigh equation (6.11) becomes singular at the cross-stream location $y^{\ast } = y^{\ast }_c$ where the flow velocity is equal to the wave velocity, $\bar {v}_x^{{{\dagger}} } = c_R^{{{\dagger}} }$. It is necessary to consider viscous effects in a ‘critical layer’ in the vicinity of the location $y^{\ast } = y^{\ast }_c$ to resolve this singularity. Due to the condition $\bar {v}_{xL}^{{{\dagger}} } < c_R^{{{\dagger}} } < \bar {v}_{xU}^{{{\dagger}} }$ (seventh row in table 2), there always exists an internal critical layer where viscous stress is significant, and so the flow is destabilised by viscosity even in the limit of high Reynolds number, as first suggested by Reynolds (1883).

The two solutions near the regular singular point can be determined using a Forbenius expansion in $(y^{\ast } - y^{\ast }_c)$:

(6.24)\begin{equation} \left.\begin{array}{c@{}} \tilde{v}^{{{\dagger}}}_{y1}(y^{{\ast}}) = (y^{{\ast}} - y^{{\ast}}_c) P_1(y^{{\ast}} - y^{{\ast}}_c), \\ \tilde{v}^{{{\dagger}}}_{y2}(y^{{\ast}}) = P_2(y^{{\ast}} - y^{{\ast}}_c) + (y^{{\ast}} - y^{{\ast}}_c) \log{(y^{{\ast}} - y^{{\ast}}_c)} P_1(y^{{\ast}} - y^{{\ast}}_c). \end{array}\right\} \end{equation}

The solution $\tilde {v}^{{{\dagger}} }_{y2}$ has a term proportional to $\log {(y^{\ast } - y^{\ast }_c)}$, which is not only singular at $y^{\ast } = y^{\ast }_c$, but is also multi-valued near the critical point. If we consider $\log {(y^{\ast } - y^{\ast }_c)} = \log {(|y^{\ast } - y^{\ast }_c|)}$ for $y^{\ast } > y^{\ast }_c$, there are two possibilities $\log {(y^{\ast } - y^{\ast }_c)} = \log {(|y^{\ast } - y^{\ast }_c|)} \pm \imath {\rm \pi}$ for $y^{\ast } < y^{\ast }_c$. To determine the correct solution, it is necessary to include the viscous terms and solve the complete Orr–Sommerfeld equation in a thin region of thickness ${Re}^{-1/3}$ around the critical point, where viscous and inertial effects are comparable, as discussed in chapter 4 of Drazin & Reid (1981). The matched asymptotic analysis was first carried out in the classic works of Tollmien (1929) and Lin (1945), and it turns out that the correct solution is $\log {(y^{\ast } - y^{\ast }_c)} = \log {(|y^{\ast } - y^{\ast }_c|)} - \imath {\rm \pi}$ for $y^{\ast } < y^{\ast }_c$, if $(\textrm {d} \bar {v}^{\ast }_x / \textrm {d} y^{\ast }) > 0$ at $y^{\ast } = y^{\ast }_c$.

Though the velocity eigenfunctions for these inviscid modes have the same mode structure as the Tollmien–Schlichting modes in a rigid-walled channel, these are not, in general, continuations of the Tollmien–Schlichting modes. For a rigid-walled channel, the only velocity scale for the Tollmien–Schlichting mode is the flow velocity, and wave speed of the Tollmien–Schlichting mode decreases to zero as the Reynolds number is increased. In contrast, there is an additional velocity scale, the shear wave velocity of the wall, in the flow through soft-walled conduits. Thus there are additional inviscid modes with a critical layer in a soft-walled conduit, which are not continuations of the rigid-wall Tollmien–Schlichting mode, but whose wave speed scales with the solid shear wave velocity. These modes do become unstable in the high-Reynolds-number limit. Thus, in addition to the modification of the rigid-wall Tollmien–Schlichting modes, there are additional inviscid modes with or without a critical layer which could become unstable in a compliant conduit.

6.3. Wall layer

At high Reynolds number, there are solutions of the linear stability equations where the viscous stress is confined to a region of thickness ${Re}^{-1/3}$ at the wall. When the distance from the wall is small compared to the characteristic length, the fluid mean velocity can be approximated as a linear profile:

(6.25)\begin{equation} \bar{v}_x^{{{\dagger}}} = \dot{\gamma}_w^{{{\dagger}}} y^{{\ast}}, \end{equation}

where $\dot {\gamma }_w^{{{\dagger}} }$ is the strain rate at the wall for the mean flow scaled by $(V_f/h_f)$. The scaled wall coordinate is defined as $y_w = \epsilon y^{\ast }$, where $\epsilon \ll 1$. This is substituted into the mass and streamwise momentum equations (2.7) and (2.8):

(6.26)$$\begin{gather} \frac{1}{\epsilon} \frac{\textrm{d} \tilde{v}_y^{{{\dagger}}}}{\textrm{d} y_w} + \imath k^{{\ast}}_x \tilde{v}_x^{{{\dagger}}} = 0, \end{gather}$$

(6.27)$$\begin{gather}\imath k^{{\ast}}_x (\epsilon \dot{\gamma}_w^{{{\dagger}}} y_w - c^{{{\dagger}}}) \tilde{v}_x^{{{\dagger}}} + \tilde{v}_y^{{{\dagger}}} \frac{\textrm{d} \bar{v}_x^{{{\dagger}}}}{\textrm{d} y^{{\ast}}} ={-} \imath k^{{\ast}}_x \tilde{p}^{{{\dagger}}} + \frac{1}{{Re} \epsilon^{2}} \left( \frac{\textrm{d}^{2}}{\textrm{d} y_w^{2}} - \epsilon^{2} k^{{\ast} 2}_x\right) \tilde{v}_x^{{{\dagger}}}. \end{gather}$$

The wall-layer velocities are defined as $\tilde {v}_{wx} = \tilde {v}_x^{{{\dagger}} }$ and $\tilde {v}_{wy} = \epsilon \tilde {v}_y^{{{\dagger}} }$ for a balance between the terms in the mass conservation equation (6.26). In the limit of high Reynolds number, there is a balance between the inertial and viscous terms in the $x$ momentum equation for $\epsilon = {Re}^{-1/3}$, $c^{{{\dagger}} } = {Re}^{-1/3} c_w$ and $\tilde {p}^{{{\dagger}} } = {Re}^{-1/3} \tilde {p}_{fw}$. Neglecting terms $O({Re}^{-1/3})$ and smaller, the scaled $x$ momentum equation is

(6.28)\begin{equation} \imath k^{{\ast}}_x (\dot{\gamma}_w^{{{\dagger}}} y_w - c_w) \tilde{v}_{wx} + \dot{\gamma}_w^{{{\dagger}}} \tilde{v}_{wy} ={-} \imath k^{{\ast}}_x \tilde{p}_{fw} + \frac{\textrm{d}^{2} \tilde{v}_{wx}}{\textrm{d} y_w^{2}}. \end{equation}

As is usual for boundary layer flows, the leading-order $y$ momentum equation states that the pressure gradient across the wall layer is zero:

(6.29)\begin{equation} \frac{\textrm{d} \tilde{p}_{fw}}{\textrm{d} y_w} = 0. \end{equation}

Equation (6.28) is differentiated with respect to $y_w$, and the mass conservation equation and zero pressure gradient condition (equations (6.26) and (6.29)) are used to obtain one equation for $\tilde {v}_{wx}$:

(6.30)\begin{equation} \imath k^{{\ast}}_x (\dot{\gamma}_w^{{{\dagger}}} y_w - c_w) \frac{\textrm{d} \tilde{v}_{wx}}{\textrm{d} y_w} - \frac{\textrm{d}^{3} \tilde{v}_{wx}}{\textrm{d} y_w^{3}} = 0. \end{equation}

The solution of (6.30) for $\tilde {v}_{wx}$ is the generalised Airy function, ${Ai}((\imath k^{\ast }_x \dot {\gamma }_w^{{{\dagger}} })^{1/3}(y_w - c_w/\dot {\gamma }_w^{{{\dagger}} }), 1)$. There is a second linearly independent solution, ${Bi}((\imath k^{\ast }_x \dot {\gamma }_w^{{{\dagger}} })^{1/3}(y_w - c_w/\dot {\gamma }_w^{{{\dagger}} }), 1)$, which diverges for $y_w \gg 1$, and the coefficient of this is set equal to zero so that the solution decreases to zero for $y_w \gg 1$. For a rigid surface, the wave velocity $c_w$ is determined from the no-slip condition $\tilde {v}_{wx} = 0$ at $y_w = 0$, and it is found (Corcos & Sellars 1959; Gill 1965) that these modes are always stable.

The Airy function solution of (6.30) is not complete, because it does not satisfy the zero normal velocity boundary condition at the surface. The inviscid part of the solution, $(\tilde {v}_{ox}, \tilde {v}_{oy}, \tilde {p}_{fo})$, is obtained by solving the Rayleigh equation (6.11) with appropriate boundary conditions at the upper wall of the channel or the centre of the tube. This solution is designated the outer solution with subscript $o$. From the zero normal velocity boundary condition, it is evident that $\tilde {v}_{oy} \sim \tilde {v}_{wy} \sim {Re}^{-1/3} \tilde {v}_{wx}$. From the mass conservation equation, the streamwise and cross-stream velocity perturbations in the outer layer are of equal magnitude, and the pressure scales as $\rho V_f \tilde {v}_{ox} \sim {Re}^{-1/3} \rho _f V_f \tilde {v}_{wx}$, as shown in table 3.

Table 3. The magnitude of the velocity and pressure in the wall layer and outer flow and the displacement and pressure in the solid (right-hand column) relative to the tangential velocity in the wall layer for the wall mode instability.

The solid dynamics is governed by (3.39)–(3.42). The characteristic length scale for the solid displacement field is the solid thickness, $h_s$, and the displacement fluctuations in the streamwise and cross-stream directions are of equal magnitude in an expansion in ${Re}^{-1/3}$. The viscous terms are neglected in an expansion in ${Re}^{-1/3}$, since they are $O({Re}^{-1})$ smaller than the inertial and elastic terms. For the stability calculation, it is necessary to calculate the leading order and the $O({Re}^{-1/3})$ correction for the displacement field. The displacement field is expressed as $\tilde {u}^{\ast }_i = \tilde {u}^{(0)}_i + {Re}^{-1/3} \tilde {u}^{(1)}_i$. Similarly, the wave velocity is expanded in a series $c_w = c_w^{(0)} + {Re}^{-1/3} c_w^{(1)}$.

The relative magnitudes of the fluid velocity and solid displacement fields are determined by the matching conditions at the interface. The magnitudes of perturbations to the outer flow velocity and solid pressure and the pressure and displacement, summarised in table 3, are briefly explained here. The elastic stress in the solid and the viscous stress in the fluid are balanced for $G \sim {Re}^{-1/3} \rho _f V_f^{2}$ or ${Re} \sim \varSigma ^{3/5}$. The asymptotic analysis of Kumaran (1998) revealed that the wall modes are always stable in this regime, but there was an instability when one of these modes was numerically continued to lower values of elasticity. A modified asymptotic analysis where ${Re} \sim \varSigma ^{3/4}$ revealed that the flow could be unstable (Shankar & Kumaran 2001a, 2002). The magnitudes of the perturbation quantities are shown in table 3, based on the assumption that the magnitude of the velocity in the wall layer is $V_f \tilde {v}_{wx}$. Here, the scaling between $\varSigma$ and ${Re}$ is determined by balancing the pressure in the solid with the pressure in the outer flow in the fluid, ${Re}^{-1/3} \rho _f V_f^{2} \tilde {v}_{wx} \sim G {Re}^{1/3} \tilde {v}_{wx}$, which reduced to ${Re} \sim \varSigma ^{3/4}$.

Using the scalings in table 3, the leading order and the first correction to the normal velocity condition at the interface (3.43a,b) are

(6.31a,b)\begin{equation} 0 = \tilde{u}^{(0)}_y \quad \text{and} \quad \tilde{v}_{wy} + \tilde{v}_{oy} ={-} \imath k^{{\ast}}_x (c_w^{(0)} \tilde{u}^{(1)}_y + c_w^{(1)} \tilde{u}^{(0)}_y), \end{equation}

and the tangential velocity conditions at the interface (3.43a,b) are

(6.32a,b)\begin{equation} \tilde{u}^{(0)}_y = 0\quad \textrm{and} \quad \tilde{v}_{wx} + \dot{\gamma}_w^{{{\dagger}}} \tilde{u}^{(1)}_y ={-} \imath k^{{\ast}}_x c_w^{(0)} \tilde{u}^{(0)}_y. \end{equation}

The leading order and first correction to the normal stress balance, (3.44), are

(6.33a,b)\begin{equation} - \tilde{p}^{(0)}_f ={-} \tilde{p}^{(0)}_s + 2 \frac{\textrm{d} \tilde{u}^{(0)}_y}{\textrm{d} Y^{{\ast}}}\quad \text{and} \quad - \tilde{p}^{(1)}_f ={-} \tilde{p}^{(1)}_s + 2 \frac{\textrm{d} \tilde{u}^{(1)}_y}{\textrm{d} Y^{{\ast}}}. \end{equation}

The leading order and first correction to the shear stress balance, (3.45), are

(6.34a,b)\begin{equation} 0 = \frac{\textrm{d} \tilde{u}^{(0)}_x}{\textrm{d} Y^{{\ast}}} + \imath k^{{\ast}}_x \tilde{u}^{(0)}_y\quad \text{and} \quad \frac{\textrm{d} \tilde{v}_{wx}}{\textrm{d} y^{{\ast}}} = \frac{\textrm{d} \tilde{u}^{(1)}_x}{\textrm{d} Y^{{\ast}}} + \imath k^{{\ast}}_x \tilde{u}^{(1)}_y. \end{equation}

The reasons for the above balances are as follows. The relation ${Re} \sim \varSigma ^{3/4}$ was determined from the condition that the normal stress in the solid and the pressure due to the outer flow in the fluid are of equal magnitude. Using this scaling, it is easily verified from table 3 that the shear stress in the wall layer, which is the dominant contribution to the fluid stress, is $O({Re}^{-1/3})$ smaller than that in the solid. Therefore, the shear stress in the solid is set equal to zero in the leading approximation in (6.34a,b), and the first correction to the shear stress in the solid is set equal to the shear stress in the wall layer in the fluid. From table 3, it is evident that the normal velocity in the fluid $\tilde {v}_y$ is ${Re}^{-1/3}$ smaller than that in the solid, $k_x c \tilde {u}_y$. Therefore, the normal displacement of the solid is set equal to zero in the leading approximation in (6.31a,b), and the first correction to the displacement field is balanced by the leading-order normal displacement in the fluid. The tangential velocity boundary condition, (3.43a,b), is interesting, because the largest term turns out the be the change in the fluid mean velocity due to the surface displacement, $\dot {\gamma }_w \tilde {u}_y$. This is $O({Re}^{1/3})$ larger the $x$ velocity of the surface, $c \tilde {u}_y$, because the wave speed $c$ is ${Re}^{-1/3}$ smaller than $V_f$. Similarly, $\dot {\gamma }_w \tilde {u}_y$ is also $O({Re}^{1/3})$ larger than the fluid velocity perturbation at the interface.

Both the normal and tangential velocity conditions prescribe $\tilde {u}^{(0)}_y = 0$ at the interface, and the leading-order tangential stress is zero from (6.34a,b). Thus, the gel dynamics in the leading approximation is that of an elastic solid with zero normal displacement and zero tangential stress prescribed at the interface. The perturbations are neutrally stable in this leading approximation, since there is no dissipation included. Therefore, $c_w^{(0)}$ is real, and there are multiple solutions for $c_w^{(0)}$ which are like the fundamental modes and harmonics of a solid layer. The natural frequencies of these modes of oscillation are $- k^{\ast }_x c_w^{(0)}$. In order to determine the stability of the system, it is necessary to calculate the first correction to the wave speed $c_w^{(1)}$. The results of the analysis show that the flow becomes unstable when the Reynolds number exceeds a critical value ${Re}_c$, which is proportional to $\varSigma ^{3/4}$ (Kumaran 1998; Shankar & Kumaran 2001a, 2002). The mechanism of destabilisation is the transfer of energy from the mean flow to the fluctuations due to the shear work done at the interface.

6.4. Numerical linear stability analyses

The flow in a channel bounded by spring-backed walls, with configuration shown in figure 1(b), was considered by Rotenberry (1992), Rotenberry & Saffman (1990), Davies & Carpenter (1997a,b) and Larose & Grotberg (1997). In these studies, with the exception of Larose & Grotberg (1997), the normal wall motion was described by a constitutive relation similar to (3.1), while the tangential velocity was zero at the wall. The effect of the compliant wall on the Tollmien–Schlichting instability was considered in Rotenberry (1992), Rotenberry & Saffman (1990) and Davies & Carpenter (1997b), who found that wall compliance has a stabilising effect which leads to transition delay. This conclusion should be interpreted with care, because it is known that the transition in a rigid channel is highly subcritical; even though linear stability analysis predicts that the transition Reynolds number is about 5772, the transition is observed in experiments at a much lower Reynolds number of about 1200 due to the highly subcritical nature of the transition. An interesting result of Rotenberry & Saffman (1990) is that the bifurcations transform from subcritical to supercritical as the wall compliance is increased. This result suggests that the linear stability analysis may be more accurate for a compliant tube than for a rigid tube; however, there does not appear to be further work on this issue.

The effect of wall compliance on the FISI was considered by Davies & Carpenter (1997a). They found two types of instabilities, travelling wave flutter and static divergence, which are also found in the boundary layer flow past a compliant wall (Carpenter & Garrad 1985, 1986). The travelling wave flutter is associated with an internal critical layer within the fluid which generated a phase shift between the normal velocity and pressure at the wall, in a manner similar to that for the Tollmien–Schlichting instability. It was shown by Huang (1998) that the flutter instability is adequately captured by an inviscid analysis, if the phase shift in the normal velocity at the critical layer is included in the analysis. The asymptotic analysis of the viscous critical layer, carried out by Gajjar & Sibanda (1996), showed that while wall compliance was stabilising, dissipation in the wall could have a destabilising effect on the travelling wave flutter. The asymptotic basis of the static divergence has not been fully understood. Tangential wall motion was included in the analysis of Larose & Grotberg (1997), using a model similar to (3.2), who considered the channel configuration shown in figure 1(b) as a model for a collapsed lung airway. The Reynolds number was lower than the critical Reynolds number for the Tollmien–Schlichting mode. In addition to the travelling wave flutter, Larose & Grotberg (1997) also reported a long-wave instability; the latter is present only when wall tangential motion is included. They also found that the travelling wave flutter and the long-wave instability could appear simultaneously, resulting in a codimension 2 bifurcation.

The stability of the pressure-driven flow in a channel with one wall made of a massless tensioned membrane was studied by Stewart et al. (2010). In this case, the membrane response was represented by an equation similar to (3.1) with zero inertia and spring constant, and the tangential velocity of the fluid at the wall was set equal to zero. Several modes of instability were identified, including the Tollmien–Schlichting mode and the travelling wave flutter. The authors found that the downstream travelling Tollmien–Schlichting mode has an asymptotic structure which is different from that in a rigid channel. Stewart et al. (2009) considered a situation where only a finite section of one of the channel walls was made of a massless tensioned membrane, and found global instabilities with characteristic mode shapes, even in situations where a local stability analysis indicated that there is no absolute instability. In contrast, the study of Kumaran & Srivatsan (1998) considered the configuration shown in figure 2(a) where the flow is driven by moving walls. Here, analytic continuation was used to extend the low-Reynolds-number viscous instability to higher Reynolds numbers. It was found that as the Reynolds number is increased, the system first becomes unstable at a transition Reynolds number, and then becomes stable again as the Reynolds number is increased. The transition Reynolds number at the lower stability boundary increases proportional to the square of the parameter $\varSigma _T = (\rho _f T h_f / \mu _f^{2})$ at low Reynolds number, where T is the surface tension, while the upper boundary is a constant.

An issue that needs to be carefully considered in the flows bounded by spring-backed plates and membranes is the pressure gradient and the shear stress acting at the wall, which results in wall deformation even at steady state. In spring-backed plate walls, the pressure gradient will result in the narrowing of the channel/tube with downstream distance, which is typically not incorporated in the stability analysis. It is necessary to assume that this narrowing is much smaller than the channel width. The shear stress exerted on the wall results in a variation in the tension with downstream distance. It is necessary to assume that the pretension is much larger than the variation in tension due to the shear stresses. These difficulties do not exist in the case of an incompressible viscoelastic wall, where it is shown in § 3.3.3 that the base state is steady and fully developed in which the displacement field is independent of the streamwise coordinate.

The Couette flow past a viscoelastic surface, shown in figure 3(a), does not contain an internal critical layer; here, the only possible mode of destabilisation is the wall mode. The high-Reynolds-number flow was studied by Chokshi & Kumaran (2008b), who identified the presence of the wall mode instability. That study included a weakly nonlinear analysis, where the domains of subcritical and supercritical bifurcations were identified. It was found that the instability is supercritical over most of parameter space, though there are small domains where the instability is subcritical. This indicates that there is a continuous transition in the flow past a compliant surface, in contrast to the flow through rigid conduits where the transition is supercritical. Pierce (1992) conducted the initial study of the Poiseuille flow in a channel bounded by viscoelastic walls of finite thickness using variational principles. The linear wall model ((3.11) without the last term on the right-hand side) was used, and the coupled equations for the wall displacement and fluid velocity were solved. A Ginzburg–Landau equation was solved for the amplitude of the interface perturbations. The boundary conditions do not appear to be complete, since only the normal velocity and stress conditions are balanced at the interface. This was perhaps the first theoretical study to suggest that the transition Reynolds number for a channel with a viscoelastic wall is lower than that for a rigid channel, and the transition Reynolds number decreases as the wall thickness increases. However, only the modification of the rigid-wall instability due to wall compliance was considered.

The stability of the parabolic flow in a tube with an annular viscoelastic wall was studied by Kumaran (1996) in the high-Reynolds-number limit for ‘regular inviscid’ modes which do not contain a critical layer. The linear viscoelastic model was used for the wall dynamics. The system was found to be neutrally stable in the inviscid approximation, but the inclusion of viscous effects stabilised the perturbations. Shankar & Kumaran (1999) carried out a similar study for a non-parabolic flow in the entrance section of a tube or in a converging tube. In all cases, the ‘regular inviscid’ modes, which do not contain a critical layer, are neutrally stable in the inviscid limit, but become stable due to viscous dissipation. However, ‘singular inviscid’ modes, which do contain a critical layer, could become unstable if the Reynolds number is increased beyond a transition value. In Shankar & Kumaran (2000), it was found that the parabolic flow is always stable for axisymmetric modes, but it does become unstable to non-axisymmetric modes due to the presence of a viscous critical layer within the flow. Non-axisymmetric modes with meridional wavenumber $n=1$ were the most unstable.

The studies of Kumaran (1996) and Shankar & Kumaran (1999, 2000) were carried out using the linear elastic model for a solid with a linear strain measure ((3.11) without the last nonlinear term on the right-hand side), which is not material frame invariant. The stability for channel and pipe flows for the neo-Hookean model was analysed by Gaurav (2009) and Gaurav & Shankar (2010). Both the singular inviscid modes with scaling ${Re}_t \propto \varSigma ^{1/2}$ and the wall modes with scaling ${Re}_t \propto \varSigma ^{3/4}$ were identified in these studies. The results for the linear and the neo-Hookean model were qualitatively similar for the most unstable mode, but there were numerical differences in the transition Reynolds number. The effect of different wall models on the stability has been compared by Patne et al. (2017) and Patne & Shankar (2019a), who also find that the high-Reynolds-number inviscid and wall modes exist for the neo-Hookean and the Mooney–Rivlin models. Thus, the high-Reynolds-number instabilities are not sensitive to the wall model, though there are quantitative differences in the transition Reynolds number for different wall models. The reason for this is discussed at the end of § 3.3, where it is shown that the strain in the base state in the solid decreases as an inverse power of the Reynolds number: ${Re}^{-1}$ for the inviscid modes and ${Re}^{-1/3}$ for the wall modes. If the terms multiplying the strain in the base state are neglected, the linear model for the viscoelastic wall is obtained.

7. Non-Newtonian fluids

There have been a few studies which have considered viscoelastic fluids (Chokshi & Kumaran 2007; Chokshi et al. 2015), where the Oldroyd B model is usually used, and the stress tensor is written as the sum of a viscous part $\sigma ^{v}$ and an elastic part $\sigma ^{p}$. The elastic stress is related to the second-order polymer ‘conformation tensor’ ${\boldsymbol c}$, which represents the dyadic product of the end-to-end vector of the polymer chains. The equation for the conformation tensor is

(7.1)\begin{equation} \frac{{\mathcal{D}} {\boldsymbol c}}{{\mathcal{D}} t} ={-} \frac{{\boldsymbol c} - {\boldsymbol c}^{eq}}{\lambda}, \end{equation}

where $\lambda$ is the polymer relaxation time and the equilibrium conformation tensor ${\boldsymbol c}^{eq} = (k_B T / H) {\boldsymbol I}$, where $k_B$ and $T$ are the Boltzmann constant and absolute temperature, $H$ is the polymer ‘spring constant’ and ${\boldsymbol I}$ is the identity tensor. The upper convected derivative $({\mathcal {D}} / {\mathcal {D}} t)$ is the rate of change of the confirmation tensor in a coordinate system moving, rotating and stretching with the fluid:

(7.2)\begin{equation} \frac{{\mathcal{D}}{\boldsymbol c}}{{\mathcal{D}} t} = \frac{\partial {\boldsymbol c}}{\partial t} + {\boldsymbol v} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol c} - {\boldsymbol c} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v} - (\boldsymbol{\nabla} {\boldsymbol v})^{T} \boldsymbol{\cdot} {\boldsymbol c}.\end{equation}

The polymer contribution to the stress tensor is then calculated from the conformation tensor:

(7.3)\begin{equation} {\boldsymbol \sigma}^{p} = \frac{\mu_p H}{\lambda k_B T} ({\boldsymbol c} - {\boldsymbol c}^{eq}), \end{equation}

where $\mu _p$ is the polymer contribution to the viscosity. The total stress is written as the sum of the polymeric and solvent stresses:

(7.4)\begin{equation} {\boldsymbol \sigma}^{f} = {\boldsymbol \sigma}^{p} + {\boldsymbol \sigma}^{v}, \end{equation}

where ${\boldsymbol \sigma }^{v}$ is the solvent stress described by Newton's law of viscosity and $\mu _f$ is the solvent viscosity. The ‘retardation’ parameter $\beta$, which is defined as $(\mu _f/\mu ) = (\mu _f/(\mu _f + \mu _p))$, is a measure of the solvent contribution to the total stress.

There are two dimensionless parameters in the Oldroyd model: the Weissenberg number, which is usually defined as ${We} = (\lambda V_f/h_f)$, and the retardation parameter $\beta$. It should be noted that the Weissenberg number is defined differently in Chokshi & Kumaran (2007) and Chokshi et al. (2015) as $(\lambda G / \mu _f)$, but the standard definition in terms of the fluid strain rate is used here. The upper convected Maxwell model is a simplification of the Oldroyd model for $\beta = 0$. The Oldroyd model in (7.1)–(7.4) does not consider spatial variations in the polymer concentration field and the flux of polymers due to concentration gradients and stresses. A more detailed model, the ‘diffusive Oldroyd’ model (Chokshi & Kumaran 2007), has been used, which couples the polymer concentration and conformation fields to the stress tensor.

Shankar & Kumar (2004) studied the stability of the Couette flow of a viscoelastic fluid past a neo-Hookean surface in the configuration shown in figure 3(a). The upper convected Maxwell model was used for the fluid, and the zero-Reynolds-number flow was considered. Earlier work by Gorodtsov & Leonov (1967) showed that there are two distinct modes for an upper convected Maxwell fluid in a channel with rigid walls which are always stable. Shankar & Kumar (2004) found that one of the Gorodtsov–Leonov modes does become unstable as the wall compliance is increased. However, viscoelasticity has a stabilising effect on the instability due to a compliant wall, and there is a critical Weissenberg number beyond which the flow is always stable.

The analysis was extended to the Oldroyd B model for a viscoelastic fluid by Chokshi & Kumaran (2007). For relatively dilute polymeric solutions with retardation parameter $\beta > 0.5$, the results were in agreement with those of Shankar & Kumar (2004) that an instability could exist only when the Weissenberg number is below a maximum value. However, for relatively concentrated solutions with $\beta < 0.5$, there could be an instability in the viscous limit for any Weissenberg number, when the scaled velocity $\varGamma = (V_f \mu _f / G h_f)$ exceeds a critical value. The weakly nonlinear analysis carried out by Chokshi & Kumaran (2007) showed that the bifurcation is subcritical when there is no dissipation in the solid, but it transforms into a supercritical bifurcation for $(\sqrt {\mu _r} / H) > 1$. Here, $\mu _r$ and $H$ are the ratio of viscosities and thicknesses of the solid and fluid.

The high-Reynolds-number wall mode instability of the Couette flow of an Oldroyd B fluid was examined in Chokshi et al. (2015). The transition Reynolds number was found to scale as ${Re}_t \propto \varSigma ^{3/4}$, which is the same as that for a Newtonian fluid, but the transition Reynolds number was found to be lower, by up to a factor of 10, in comparison with a Newtonian fluid. Thus, it was found that viscoelasticity has a strong destabilising effect on the wall mode instability. However, in a small region in parameter space for very dilute polymers with $\beta > 0.9$ and when the solid thickness is less than the fluid thickness, viscoelasticity had a stabilising effect. Thus, by a careful choice of parameters, it is possible to stabilise or destabilise the flow. Numerical calculations also indicated the presence of an inviscid instability with transition Reynolds number ${Re}_t \propto \varSigma ^{1/2}$ which could be smaller than that for the wall mode at high Reynolds number. However, the growth rate of this instability was found to be smaller than that of the wall mode instability by a factor of $10^{-3}$, suggesting that these might be difficult to observe in experiments. The decrease in the transition Reynolds number due to viscoelasticity was observed in the experiments of Srinivas & Kumaran (2017a) discussed in § 8.

The power-law and Bingham models are the other class of models used for non-Newtonian fluids. The stress has the form

(7.5)\begin{equation} {\boldsymbol \sigma} = {\boldsymbol \sigma}_y + \kappa (\boldsymbol{\nabla} {\boldsymbol v}) |\boldsymbol{\nabla} {\boldsymbol v}|^{n-1}, \end{equation}

where ${\boldsymbol \sigma }_y$ is the yield stress and $n$ is the power-law index. For the Bingham model, the power-law index is $1$ and $\kappa$ is the fluid viscosity. For the power-law model, the yield stress is zero. Patne & Shankar (2019b) studied the stability of the pressure-driven flow of a Bingham fluid, where the stress is the sum of a yield stress and a viscous stress, through a channel with compliant walls. This was different from the results for the same flow of a Newtonian fluid (Gaurav & Shankar 2010; Patne et al. 2017), where no instability was found. There is an unyielded region at the centre of the channel in a Bingham fluid which is not present for a Newtonian fluid, and the authors suggested that the unyielded region acts as a moving upper wall for the flowing region, in a manner similar to that for the Couette flow in a channel. This resulted in a finite-wavelength instability similar to that observed for the Couette flow in a channel in the viscous limit (Kumaran et al. 1994; Chokshi & Kumaran 2008a).

The stability of wall modes at high Reynolds number in the Couette flow of a power-law fluid past a compliant surface was studied by Giribabu & Shankar (2017). The authors carried out high-Reynolds-number asymptotic and numerical studies. The asymptotic analysis showed that the thickness of the viscous wall layer scales as ${Re}^{-n/(2n+1)}$ times the channel width; this reduced to the scaling ${Re}^{-1/3}$ for a Newtonian fluid with $n=1$. The transition Reynolds number scales as ${Re}_t \propto \varSigma ^{(2n+1)/2(n+1)}$. The numerical results were in agreement with the asymptotic analysis when the power-law index was greater than $0.3$, but the scaling of the transition Reynolds number was closer to that for the inviscid modes, ${Re}_t \propto \varSigma ^{1/2}$, for power-law index less than $0.3$.

8. Experimental studies

In most experimental studies, the soft walls for channels/tubes have been fabricated using two materials, polyacrylamide and polydimethylsiloxane (PDMS) gels. Polyacrylamide gels are used in gel electrophoresis for separation of biological molecules based on their electrophoretic mobility. These are transparent hydrophilic materials which consist of a three-dimensional cross-linked polymer network swollen with water. The acrylamide molecules are cross-linked by methylene bisacrylamide cross-linker, and ammonium persulphate and tetramethylene diamine act as initiator and terminator for the chain reaction. The cross-linking of the acrylamide monomers takes place at room temperature in a few hours. The cross-link density and the elasticity modulus are varied by changing the proportions of the different constituents, and shear elasticity modulus can be varied in the range 1 kPa to 0.5 MPa. The bulk elasticity modulus is 10–50 times the shear modulus. The channels/tubes are made by template-assisted fabrication, where a rod or a slide made of a hard material like glass is used as the template. A container is fabricated in the shape of the outer surface of the channel wall, and the template is held in position at the centre of the container. The pre-gelation mixture is poured into the container, and the cross-linking takes place at room temperature. After the gelation is complete, the template is withdrawn, resulting in a conduit of rectangular or circular cross-section within the gel block. The characteristic dimension of the conduit has to be 0.5 mm or more, since a feature of size of less than 0.5 mm may not retain its integrity.

A PDMS gel is a silicone-based gel which is used in microfluidics and lab-on-a-chip applications. This gel is also transparent, but it is hydrophobic and is not wetted by water. The monomer is dimethylsiloxane with dimethylvinyl termination, and the cross-linker is usually a dimethylhydrogensiloxane. A catalyst, chlorplatinic acid which is a platinum complex, and an elevated temperature of about 60–80 $^{\circ }$C are required for the cross-linking reaction. Polydimethylsiloxane is used for template-assisted fabrication, but is especially suitable for soft lithography, where a pattern is transferred from a hard substrate to a soft gel. The steps used in a typical process for fabricating microchannels is shown in figure 7. Usually, the negative of the desired channel is patterned on a silicon wafer using a photoresist such as SU-8. The photoresist of the desired thickness is spin-coated on to a silicon wafer. This is exposed to light in the patterned regions, where the photoresist is chemically modified usually by cross-linking. The non-exposed regions are washed off using a developer, resulting in a photoresist pattern on the wafer in the exposed regions. The PDMS pre-gelation mixture is poured onto the wafer, and the PDMS is cross-linked in an oven. The PDMS stamp is then peeled off the silicon wafer, containing the channel features as indentations in the PDMS block. This is then closed, usually with a glass plate or another PDMS layer, using a process called plasma bonding. Here, the PDMS is treated with oxygen plasma, which exposes silanol groups on the surface. When contacted with a glass or another PDMS layer similarly exposed, the formation of Si–O–Si bonds between the two layers results in strong adhesion and a leak-proof channel. The procedure for PDMS fabrication is more complicated than that for polyacrylamide fabrication. The minimum elasticity modulus of about 10 kPa is higher than that for polyacrylamide. The important advantage of PDMS is that the minimum feature size that can be transferred to a PDMS surface is of the order of micrometres, and so it is possible to fabricate features of very small dimensions.

Figure 7. Steps in a typical fabrication procedure for a microchannel in PDMS using soft lithography.

The pioneering experiments on the effect of wall compliance on the transition in a tube were performed by Lahav et al. (1973) and Krindel & Silberberg (1979). The tube was a cylindrical bore of diameter about 0.15 mm at the centre of a cylinder of polyacrylamide gel of length 5 cm and diameter 2.6 or 9 mm. A cylindrical cavity of diameter 2.6 or 9 mm was drilled into a rigid Perspex block, the pre-gelation reaction mixture was poured into the block, a hypodermic needle of 30 gauge was inserted into the centre and the ends were closed. After completion of the cross-linking, the needle at the centre of the tube was withdrawn, leaving a cylindrical bore at the centre. The ends of the tube were connected to the fluid inlet and outlet, and the pressures at the inlet and outlet were measured using pressure transducers. The authors measured the flow rate $Q$ as a function of the pressure drop ${\rm \Delta} p$ across the tube. When a pressure difference is applied across the ends, there is a variation in the tube diameter, and this variation was measured using imaging under a microscope. A provision was made to insert a dye stream at the centre of the tube, and examine the flow of the dye under a microscope. A schematic, not to scale, of the experimental set-up depicting a tubular bore inside a cylindrical gel block is shown in figure 8(a).

Figure 8. (a) Schematic, not to scale, of the experimental set-up of Krindel & Silberberg (1979) and (b) figure 6 of Krindel & Silberberg (1979) showing the scaled flow rate $(Q/Q_0)$ as a function of the Reynolds number. Here, $Q$ is the flow rate in the experiments and $Q_0$ is the flow rate for a laminar flow for the same conduit geometry and pressure drop. The horizontal axis in (b), $\log _{10}({Re})$, varies from $10^{1.4}$ to $10^{3.6}$. In (b), line (2) for the rigid tube exhibits a discontinuous transition and line (3) is for the gel-walled tube that has a lower flow rate in comparison to the rigid tube for the same Reynolds number. Panel (b) reprinted with permission from Krindel & Silberberg (1979).

The scaled flow rate was defined as $(Q/Q_0)$, the ratio of the actual flow rate and the flow rate expected for a Hagen–Poiseuille flow for the same pressure drop, $Q_0 = ({\rm \pi} {\rm \Delta} p / (8 \mu _f L_t \langle R^{-4} \rangle ))$, and the Reynolds number was based on the flow rate and tube radius, ${Re} = (2 \rho _f Q / {\rm \pi}\mu _f) \langle R^{-1} \rangle )$. Here, $L_t$ is the length of the tube, ${\rm \Delta} p$ is the pressure difference between the ends, $R$ is the local tube radius and $\langle R^{-n} \rangle$ is the average of $R^{-n}$ carried out over the length of the tube. The scaled flow rate $(Q/Q_0)$ is $1$ for a laminar flow with a parabolic velocity profile and is less than $1$ for a turbulent flow. An example of the result is shown in figure 8(b), where $(Q/Q_0)$ is plotted as a function of $\log {({Re})}$. The filled symbols are the results for a rigid tube and the open symbols are the results for gel-walled tubes with two different shear elastic moduli. The dashed line (1) is the scaled flow rate for the laminar flow, which differs from $1$ because the pressure has been adjusted for the kinetic energy of the flow. The authors report a transition at a Reynolds number of $700$ for the rigid tube, where the flow rate decreases discontinuously at a fixed pressure drop, as shown by line (2) in figure 8(b). For the flow through the gel-walled tube, the scaled flow rate is significantly lower than that for a rigid tube, and it decreases continuously. Based on visualisation of the dye stream introduced at the centre of the channel, the authors reported that there was a transition to turbulence at a Reynolds number of 570 when the elasticity modulus of the gel wall was 320 Pa and about 870 when the elasticity modulus of the gel wall was 1200 Pa. However, the scaled flow rate decreased below $1$ at a lower Reynolds number for both the shear moduli; the authors interpreted this to mean that the flow rate near the wall was turbulent even when the flow rate at the centre containing the dye stream was laminar.

The experiments of Krindel & Silberberg (1979) have been criticised for many reasons. The transition Reynolds number of $700$ for a rigid tube reported by them is much lower than the commonly accepted transition Reynolds number of 2100. The computational study of Yang et al. (2000) indicated that the decrease in the scaled flow rate can be explained on the basis of the tube deformation, without having to postulate a decrease in the transition Reynolds number. Nevertheless, Lahav et al. (1973) and Krindel & Silberberg (1979) were the first to suggest that the transition in a compliant tube could take place at a lower Reynolds number than that in a rigid tube, and that the transition was induced by wall oscillations.

8.1. Low-Reynolds-number transition

Experiments on the Couette flow past a compliant surface in the low-Reynolds-number limit have been carried out, but there have not been any experiments on pressure-driven flow. This is because soft deformable materials are unlikely to be able to withstand the pressure difference and the shear stress required to drive a viscous fluid through a channel/tube at low Reynolds number. For a low-Reynolds-number flow, the viscous stress in the fluid and the elastic stress in the solid are comparable for $(\rho _f V h_f / \mu _f) \ll 1$ and $(V \mu _f / G h_f) \sim 1$, where $V$ is the characteristic velocity and $h_f$ is the characteristic length scale in the fluid. This implies that the velocity-independent parameter $(\rho _f G h_f^{2} / \mu _f^{2}) \ll 1$. The minimum elasticity modulus for soft materials such as polymer gels is ${\sim }10^{3}$ Pa and the density of most liquids is approximately $10^{3}$ kg m$^{-3}$. If we consider a characteristic dimension of 1 mm, $(\rho _f G h_f^{2} / \mu _f^{2}) <1$ for very viscous fluids with viscosity $\sim$1 kg m$^{-1}$ s$^{-1}$ (about $10^{3}$ times the viscosity of water). The pressure gradient required to drive the flow scales as $(\mu _f V / h_f^{2})$, which is ${\sim }(G / h_f)$, if we assume the scaled velocity $(V \mu _f / G h_f)$ is $1$. Thus, the pressure difference across the length $L$ of the tube is $(G L/h_f)$, and the pressure at the inlet of the tube is larger, by approximately the ratio of the tube length and diameter, than the elasticity modulus of the tube.

In the linear stability analyses, it is assumed that the wall of the tube is made of an incompressible material. Real elastic materials are compressible, though their bulk modulus is much larger than their shear elasticity modulus. For example, Srinivas & Kumaran (2017b) report that the bulk modulus of polyacrylamide gel is 10–50 times the shear modulus. Therefore, for the flow in a channel/tube with length 10 times the tube diameter/channel width or greater, the gel wall cannot be considered incompressible. Moreover, the failure stress for polymer gels is 2–5 times the elasticity modulus, and so material failure is likely to result if the channel/tube length is 10 or more times the tube diameter/channel height. Due to these constraints, experiments in pressure-driven fully developed flows are difficult to accomplish at low Reynolds number. The development of the equivalent of superelastic materials that are highly deformable and have high failure stress is necessary to study the viscous instability in pressure-driven channel/tube flows.

The experiments on the Couette flow past a compliant surface have used one standard experimental set-up, which is a modification of the commercially available rheometer, shown in figure 9(a). The fluid layer is placed on the static base, and the top plate attached to a rotor is lowered on to the fluid layer. The diameter of the top plate is in the range 2–4 cm and the fluid layer thickness is in the range 0.3–1 mm. For the measurement of the viscosity of a Newtonian fluid, a constant torque is then applied on the rotor in the stress-controlled mode of measurement, and the angular velocity is measured. The stress and strain rate at the outer edge of the rheometer are calculated assuming that the flow is laminar, and the fluid viscosity is the ratio of the stress and the strain rate.

Figure 9. Schematic of (a) a rheometer showing the base, the rotor and the fluid and (b) a modified rheometer with a gel slab placed on the base, and the fluid above the gel slab.

For the experiments on the viscous instability in the flow past a gel, a slab of gel is placed on the base. The fluid layer is placed on the gel slab, and the top plate is lowered on to the fluid, as shown in figure 9(b). The thickness of the gel slab is in the range 1–5 mm and so the gel thickness is usually larger than the fluid thickness. The first experiments on the viscous flow past a gel surface were carried out by Kumaran & Muralikrishnan (2000) and Muralikrishnan & Kumaran (2002), followed by Eggert & Kumar (2004) on the dynamics after transition and Neelamegam, Giribabu & Shankar (2014) on single-layer and two-layer polymer gels. The compliant surface was a polyacrylamide gel slab with shear modulus in the range 1–4 kPa and the fluid was silicone oil with viscosity about 1 kg m$^{-1}$ s$^{-1}$. Since silicone oil is hydrophobic, it is not miscible with acrylamide, and so there is no swelling of the gel due to contact with the fluid. The fluid layer thickness was varied in the range 0.3–1 mm and the thickness of the gel slab was about 4.5 mm. A schematic of the measurement is shown in figure 10(a), where the viscosity (ratio of stress and strain rate) is plotted as a function of the strain rate. For a Newtonian fluid between two rigid plates, the blue curve shows that the viscosity remains a constant as the stress or the strain rate is increased. There is usually some variation at low strain rate due to instrument resolution limitations.

Figure 10. (a) The variation of viscosity with strain rate for a stress ramp over a rigid surface (blue line) and a gel surface (red line) and (b) the scaled critical strain rate $\varGamma = (\mu_f\dot {\gamma }_c/G)$ as a function of $(H/R)$ for gels with different thicknesses and elasticity moduli, where $\dot {\gamma }_c$ is the critical strain rate. The solid and dashed lines in (b) are the theoretical predictions for elasticity moduli $1$ and $4$ kPa, respectively. Panel (b) reprinted with permission from Kumaran & Muralikrishnan (2000).

For the flow of the same fluid over a gel slab, the viscosity reported by the rheometer software is a constant up to a critical value, after which it starts to increase dramatically, as shown by the red line in figure 9(b). If a stress ramp is applied in the stress-controlled mode, the strain rate decreases as the viscosity reported by the rheometer software increases. It the experiment is stopped after the viscosity increases to very large values, damage is observed on the surface of the gel. However, if the experiment is stopped before the viscosity increases significantly, there is no damage on the surface. If the experiment is started again and the shear stress is ramped up from zero, the dramatic increase in the viscosity is again observed at the same strain rate.

Recall that the rheometer software determines the viscosity from the stress and strain rate at the outer edge of the top plate assuming the flow is laminar. If there is a transition from a laminar flow to a more complicated flow profile with higher dissipation, it is expected that the viscosity reported by the rheometer will be higher than the material viscosity of the fluid. Thus, the sharp increase in the viscosity at a critical value of the strain rate is indicative of a transition from a laminar flow to a more complicated flow profile. This is compared with the theoretical prediction for the viscous instability of the Couette flow past a compliant surface. Here, the assumption is that the flow near the outer edge of the top plate, of diameter 2–4 cm, is similar to a Couette flow because the gap thickness, 0.3–1 mm, is much smaller than the plate radius.

The critical strain rate was found to be in good quantitative agreement with the theoretical predictions for the viscous instability in Kumaran & Muralikrishnan (2000) and Neelamegam et al. (2014). A comparison of the experimental result of Kumaran & Muralikrishnan (2000) with the theoretical prediction of Kumaran et al. (1994) is shown in figure 10(b). Here, the scaled critical strain rate $\varGamma$ is shown as a function of the ratio of the gel and fluid thickness $H$ for gels with different shear moduli in the range 1–4 kPa. There is quantitative agreement with no adjustable parameters. Neelamegam et al. (2014) also examined the flow past a two-layered gel, consisting of two layers with two different shear moduli $G_1$ and $G_2$, and with two different thicknesses. Transition was observed in the two-layer gel as well, and the strain rate was found to be in quantitative agreement with the predictions of the linear stability analysis for a viscous flow. The authors defined an effective shear modulus for the two-layered gel which was correlated to the strain rate for transition.

Though the study of Kumaran & Muralikrishnan (2000) and Muralikrishnan & Kumaran (2002) established that there is an instability at a critical strain rate, the nature of the flow after transition was not explored. Eggert & Kumar (2004) carried out a series of experiments in order to examine the nature of the flow after transition. The rheometer set-up was as shown in figure 9, and the fluid used was polypropylene oxide which is immiscible with the polyacrylamide gel. For a constant stress ramp, these experiments confirmed the transition in the viscosity–strain rate graphs of the kind shown in figure 10(a). The authors also programmed the rheometer to ramp the stress at a constant rate to a value just above the transition value, and then hold the stress constant. The important observation by the authors, shown in figure 11(a), is that above the transition stress, even if the stress is held constant, there are periodic oscillations. A steady state with a constant stress value was not observed after transition. The period of the oscillations, which was in the range 10–20 s, was much larger than the period of rotation of the rheometer top plate ($\sim$1 s), the inverse of the fluid strain rate ($\sim$4 ms) or the ratio of the viscosity and the elasticity modulus $(\mu _f/G) \sim 0.5$s. Thus, the post-transition flow was found to have a definite period much larger than the intrinsic or imposed time scales. The other important feature observed by Eggert & Kumar (2004) was the presence of hysteresis in the apparent viscosity–strain rate curves, as shown in figure 11(b). In these experiments, the stress is ramped up at a constant rate until the viscosity reaches a preset value, and then the stress is decreased at the same rate. The apparent viscosity–strain rate curves for stress ramps at different rates are shown by the different symbols. When the stress is increased at a constant rate, the apparent viscosity traverses the same curve independent of the rate at which stress is ramped up. When the stress is decreased, the trajectory of the apparent viscosity does depend on the stress ramp rate, and there was variation between different experimental runs. Another important observation was that even when the stress is decreased below its transition value, oscillations are observed in the viscosity; these are absent when the stress is increased. This indicates that the transition is subcritical and there could be bistable states close to transition. The flow after transition is complex, and is not just a non-viscometric steady flow, of the type observed in the Taylor–Couette or Rayleigh–Bénard instabilities. The flow is complex and time-dependent, and is likely sensitive to small variations in experimental conditions. However, the subcritical nature of the transition observed in Eggert & Kumar (2004) is in agreement with the predictions of the weakly nonlinear stability analysis of Shankar & Kumaran (2001b) and Chokshi & Kumaran (2008a).

Figure 11. The variation of the apparent viscosity with time when the stress is increased just above the critical stress and held constant (a), and the stress is increased and then decreased at different rates (b). Reprinted with permission from Eggert & Kumar (2004).

8.2. High-Reynolds-number transition

Transition occurs at a Reynolds number of about 2100 in a rigid tube and at about 1000 in a two-dimensional channel of infinite extent in the spanwise direction. An instability in a conduit due to compliant walls will be detected in experiments only if it occurs at a lower Reynolds number in comparison with that for a rigid conduit. This requirement places restrictions on the dimensions of the conduit and the flow velocity required for observing the transition in experiments. For the high-Reynolds-number inviscid modes, it is shown in § 6 that there is a balance between the inertial stress in the fluid and the elastic stress in the wall for $\rho _f V_f^{2}/G \sim 1$, and the velocity required to observe the predicted transition scales as $(G/\rho _f)^{1/2}$. The minimum elasticity modulus that can be attained in soft polymer gels is approximately 1 kPa, and the density of most materials is ${\sim }10^{3}$ kg m$^{-3}$, and therefore the minimum velocity required for observing the instability is 1 m s$^{-1}$. For this average flow velocity, the Reynolds number $(\rho _f V_f h_f / \mu _f)$ is below $1000$ only for $h_f <$1 mm if we consider a fluid such as water with $\rho _f \sim 10^{3}$ kg m$^{-3}$ and $\mu _f \sim 10^{-3}$ kg m$^{-1}$ s$^{-1}$. It is necessary to carry out experiments in conduits with dimensions 1 mm or less with soft walls, and generate a flow velocity of the order of 1 m s$^{-1}$. Similarly, for the wall mode instability, it is shown in § 6.3 that the transition Reynolds number scales as ${Re}_t \sim \varSigma ^{3/4}$ and it is necessary to design experiments such that $\varSigma ^{3/4} = (\rho _f G h_f^{2}/\mu _f^{2})^{3/4} < 1000$. If we again consider water flowing through a conduit with soft gel walls of elasticity modulus $10^{3}$ Pa, the characteristic length $h_f \lesssim 10^{-4}$ m, and the velocity is approximately 1 m s$^{-1}$. Therefore, the instabilities can be observed in practice only for channels and tubes of characteristic length $100\ \mathrm {\mu }$m–1 mm, and with relatively large velocity of the order of 1 m s$^{-1}$.

Wall deformation is an important issue in experiments on soft-walled conduits. In theoretical studies, a planar channel with parallel walls or a cylindrical tube is considered, the wall material is considered to be incompressible and the flow and solid displacement field are considered to be steady and fully developed. Due to this, wall deformation due to the streamwise variations in pressure is always present in experiments. While it is not possible to avoid wall deformation, the theoretical analysis could be modified to incorporate wall deformation. One procedure that has been used is to measure the dimensions for the deformed wall, and then compute the laminar velocity profile for the deformed shape of the channel or tube. This can then be used in the linear stability analysis to estimate the transition Reynolds number. In pressure-driven flows, the wall deformation does depend on the specific configuration and, more importantly, on the length of the conduit. Therefore, it is not possible to directly compare the experimental results with theoretical predictions for the steady fully developed base state in a rectangular channel or a cylindrical tube. The stability analysis has to be carried out for the laminar flow profile for the specific deformed geometry of the channel/tube.

Even in a deformed channel/tube, the shear strain in the wall is small. If we consider a flow with velocity 1 m s$^{-1}$ in a channel/tube with characteristic dimension 0.1–1 mm, the viscous stress at the wall is 1–10 Pa for a fluid such as water with viscosity $10^{-3}$ kg m$^{-1}$ s$^{-1}$. If the elasticity modulus of the wall material is 1–10 kPa, the shear strain in the material is $10^{-3} - 10^{-2}$. The small magnitude of the strain in the material has two implications. The first is that it will be very difficult to observe the surface displacement in experiments. Even if we consider a gel wall with a rather large thickness of 5 mm, the wall displacement due to the strain will be of the order of 5–50 $\mathrm {\mu }$m. The fluctuations in the wall displacement, which could be even smaller in magnitude, will be close to the optical resolution. The second implication is that the linear model for the stress–strain relationship in the solid would be a good approximation for high-Reynolds-number flows, and it may not be necessary to use more sophisticated frame-invariant models.

Though the shear strain is small, the compressional strain is not necessarily small in pressure-driven flows at high Reynolds number. The ratio of the pressure and wall shear stress scales as the ratio of the length to the exit and the channel tube/diameter. The ratio of the length and diameter could be in the range $10^{2} - 10^{3}$ in the tube flow experiments described in § 8.2.1, where the tube diameter is approximately 1 mm and the length is about 10 cm, while in the channel flow experiments discussed in § 8.2.2, the channel height is 100–160 $\mathrm {\mu }$m and the channel length is 3 cm. In these cases, the pressure at the entrance is $10^{2} - 10^{3}$ times the wall shear stress. Even if the compression modulus is $10$ times the shear modulus, the compressional strain is larger than the shear strain by a factor $10 - 10^{2}$. Using the estimate $10^{-3} - 10^{-2}$ for the shear strain calculated in the preceding paragraph, the compressional strain could be in the range $10^{-2} - 1$. Thus, compressional deformation perpendicular to the wall is expected to be much larger than the shear deformation along the wall. Though the compressional deformation is large, the slope of the wall $\alpha$ is still small in the experiments. The displacement perpendicular to the surface is the product of the compressional strain and the wall thickness, and the slope of the wall is the ratio of the displacement and the length of the channel/tube. Thus, the wall slope is compressional strain times the ratio of the wall thickness and length of the tube/channel, which is small in the experiments conducted so far.

The other important parameter is the characteristic frequency for the wall dynamics. For a flow dominated by inertia, the frequency scales as $(G / \rho _f h_f^{2})^{1/2}$; for a viscous-dominated flow, the frequency scales as $(G / \mu _f)$. Considering the lower bound of $G = 1$ kPa and the viscosity and density of water, the frequency for a flow dominated by viscosity is $10^{6}$ s$^{-1}$, while that for an inertia-dominated flow is $10^{3}$ s$^{-1}$ if we consider a channel/tube of characteristic dimension 1 mm. Wall motion with amplitude of the order of a few micrometres and frequency in the range $10^{3} - 10^{6}$ Hz is experimentally very difficult to detect.

Another issue of importance is the flow conditioning and accurate flow characterisation at the entrance to the channel/tube. In experimental studies on transition in rigid tubes/channels, great care is taken to ensure that the flow at the entrance is laminar and that there are no disturbances. The development length for a fully developed flow is ${\sim }0.06 {Re}$ times the tube diameter/channel height, and in experiments such as those of Patel & Head (1969), the entrance length is hundreds of times the pipe diameter. The requirements on the uniformity of the surface of the pipe/channel are also stringent – in the experiments of Darbyshire & Mullin (1995), the variations in the diameter of a 20 mm pipe were controlled at ${\pm }0.01$ mm. Such careful control is not possible in millimetre-sized tubes and channels made with flexible materials, but it is essential to ensure that there are no disturbances at the inlet to the soft conduit. For this purpose, there is usually a ‘development section’ upstream of the test section which is made of hard non-deformable material. This development section has to be of sufficient length that there is a fully developed flow at the entrance to the test section, and the disturbances introduced by the upstream flow manifold are damped out. It is also essential for there to be seamless bonding between the development and test sections, so that irregularities at the joint do not generate disturbances. In Verma & Kumaran (2012, 2013) and Neelamegam & Shankar (2015), for example, the development and test sections are made of the same gel, but with different catalyst concentrations, so that the gel in the development section is rigid with a shear modulus of about 0.5 MPa, while that in the test section has a shear modulus in the range 17–55 kPa.

8.2.1. Tube flow

Verma & Kumaran (2012), Neelamegam & Shankar (2015) and Chandra, Shankar & Das (2019) carried out experiments on soft-walled tubes fabricated using the template-assisted method with slender glass rods as the template. The walls of the tubes were made of PDMS having shear modulus in the range 17.5–86.4 kPa in all the studies, and the tube diameter was 0.8–1.2 mm in the case of Verma & Kumaran (2012), about 1.65 mm in Neelamegam & Shankar (2015) and about 0.4 mm in Chandra et al. (2019). Upstream of the test section was a development section of the same cross-section and length made of ‘hard’ PDMS with shear modulus 0.5 MPa. Transition was detected using different methods. Verma & Kumaran (2012) and Neelamegam & Shankar (2015) detected transition using the dye-stream method, where a stream of dye was injected at the inlet of the development section, and transition was detected from the break-up of the dye stream. The friction factor method was also used by Verma & Kumaran (2012) and Neelamegam & Shankar (2015), where the the pressure drop across the test section was measured with a pressure transducer. The Fanning friction factor $f$ for pipe flow is expressed in terms of the flow rate as

(8.1)\begin{equation} f = \frac{{\rm \pi}^{2} d^{5}}{32 Q^{2} \rho_f} \frac{{\rm \Delta} p}{L},\end{equation}

where $d$ is the tube diameter, $Q$ is the flow rate, ${\rm \Delta} p$ is the pressure difference across the test section and $L$ is the length of the test section. The Reynolds number based on the average velocity and tube diameter, when expressed in terms of the flow rate, is ${Re} = (4 \rho _f Q/{\rm \pi} d \mu _f)$. The friction factor is $(16/{Re})$ for a laminar flow, and transition is detected when the friction factor departs from $(16/{Re})$. For a tube that deforms due to the pressure gradient, it is appropriate to use $\langle d^{5} \rangle$ instead of $d^{5}$ in (8.1), and $\langle d \rangle$ instead of $d$ in the definition of the Reynolds number, where $\langle d \rangle$ is the diameter averaged over the streamwise coordinate. Verma & Kumaran (2012) also attempted to detect wall oscillations by directing a He–Ne laser at the interface between the fluid and the wall, and measuring the intensity of the scattered light as a function of time. The root mean square of the temporal fluctuations in the scattered light intensity was used to detect the onset of wall oscillations. Chandra et al. (2019) measured the ratio of the maximum and mean velocity of the tube, which is $2$ for a laminar flow and decreases below $2$ for a turbulent flow. Chandra et al. (2019) also measured the ratio of the root mean square of the velocity fluctuations at the onset of turbulence.

These studies reached the same broad conclusions: there is a transition at a Reynolds number as low as 500 for the tubes with the lowest shear modulus in the studies of Verma & Kumaran (2012) and Neelamegam & Shankar (2015) and as low as $250$ in the study of Chandra et al. (2019) with a smaller tube diameter. However, there were differences in the detail which indicate that the origins of the instability may not be the same. An example of the dye-stream visualisation of the flow through a tube of diameter 1.2 mm in a polymer gel of shear modulus 17 kPa from the study of Verma & Kumaran (2012) is shown in figure 12. The images of the flow of the dye stream in the tube are shown at different downstream locations – the location of the joint between the development and test section is $L=0$, the development section is on the left-hand side and the test section is on the right-hand side. It is observed that even at a Reynolds number of 1000, there is break-up of the dye stream in the downstream section of the tube, though the dye stream near the entrance of the test section is intact. As the Reynolds number increases, there is upstream progression of the break-up of the dye-stream. This is in contrast to the flow in a rigid tube, where the transition occurs at a Reynolds number of about 2100. Figure 12 also shows that there is substantial deformation of the tube in the test section. There is a significant expansion just downstream of the start of the test section, and then a decrease in the tube diameter with downstream distance. The diameter of the tube was measured at different downstream locations in the experiments, and these were used to determine the length-averaged functions of the diameter used in the friction factor (equation (8.1)) and Reynolds number.

Figure 12. Images at different downstream locations of the flow in a gel-walled tube of diameter 1.2 mm in the undeformed state made of PDMS of shear modulus 17.5 kPa in the test section (right of $L=0$) and 83 kPa in the development section (left of $L=0$). (a) Reynolds number of $1000$ and (b) Reynolds number of 1030. Reprinted with permission from Verma & Kumaran (2012).

The early transition is also observed in figure 13, where the friction factor is shown as a function of the Reynolds number. The black line is the result for a rigid tube, in which both the development and test sections have elasticity modulus 83 kPa. Here, the friction factor departs from the $(16/{Re})$ line at the rigid-tube transition Reynolds number of 2100. As the elasticity modulus is decreased, the transition Reynolds number decreases monotonically, and it has a minimum of about 900 for a tube of diameter 1.2 mm with a wall of modulus 17 kPa and about 500 for a tube of diameter 0.8 mm with a wall of modulus 17 kPa. The nature of the friction factor curves is also qualitatively different for the transition in a gel-walled tube. Whereas there is a discontinuity in the friction factor for the flow through a rigid tube at a Reynolds number of 2100, the transition is continuous for a gel-walled tube. Moreover, the friction factor exhibits an increase with Reynolds number for a gel-walled tube, in contrast to the monotonic decrease in the friction factor for a rigid tube.

Figure 13. The friction factor as a function of the Reynolds number for the flow in a tube of diameter 1.2 mm made of polyacrylamide gel with shear modulus of 0.5 MPa (black), 86 kPa (brown), 38 kPa (green), 25 kPa (blue) and 17.5 kPa (red). Redrawn, with permission, from Verma & Kumaran (2012).

In order to rule out the effect of the change in shape of the tube on the transition, Verma & Kumaran (2012) also carried out experiments using a rigid gel-walled tube, made of hard gel with shear modulus 0.5 MPa, in which the diameter varies with downstream distance in the same manner as the gel-walled tube with shear modulus 17.5 kPa. The transition in this tube of non-uniform diameter was found to be at a much higher Reynolds number than that in a gel-walled tube, indicating that the transition is due to a dynamical coupling between the flow and the wall, and not due to wall deformation alone.

In the dye-stream visualisation experiments, the amplitude and the frequency of the oscillations were measured near transition. The amplitude was found to increase continuously, in contrast to the discontinuous change found in a rigid tube (Reynolds 1883), and the frequency was in the range ($0.01 - 0.1) \times (G / \rho _f d^{2})^{1/2}$ in the experiments of Verma & Kumaran (2012) and Neelamegam & Shankar (2015). All of these results suggest that the transition in a gel-walled tube is qualitatively different from that in a rigid tube, and that the transition is due to a dynamical coupling between the flow and the wall dynamics.

The transition Reynolds number is shown as a function of the parameter $\varSigma$ in figure 14. Verma & Kumaran (2012) found that the scaling law ${Re}_t \propto \varSigma ^{5/8}$ provided the best fit for the transition Reynolds number. This is between the power laws ${Re}_t \propto \varSigma ^{1/2}$ for the inviscid modes and ${Re}_t \propto \varSigma ^{3/4}$ for the wall modes. This was explained by considering the flow in the downstream converging section, where the flow is more plug-like in comparison to the laminar flow. Due to this, the strain rate at the wall is higher than that for a laminar flow by a factor $({Re} \alpha )^{1/2}$, where $\alpha$ is the slope of the wall of the tube. The slope was estimated as $\alpha \sim ({Re}/\varSigma )$ in Verma & Kumaran (2012), and this provided the scaling law ${Re}_t \propto \varSigma ^{5/8}$ for the wall mode instability. Neelamegam & Shankar (2015) found that a scaling law ${Re}_t \propto \varSigma ^{3/2}$ provided the best fit for the transition Reynolds number, as shown in figure 14. This could not be satisfactorily explained based on the slope of the wall in the converging or diverging sections. The range of $\varSigma$ for the experiments of Neelamegam & Shankar (2015) extends only over a factor of $2$ on a logarithmic scale, which may be too small to reliably infer scaling laws. The experiments of Chandra et al. (2019) in a tube of diameter $400\ \mathrm {\mu }$m at lower Reynolds number showed a scaling law ${Re}_t \propto \varSigma ^{0.65}$, which is quite close to ${Re}_t \propto \varSigma ^{5/8}$ proposed by Verma & Kumaran (2012).

Figure 14. The transition Reynolds number as a function of the parameter $\varSigma$ from Verma & Kumaran (2012) (black), Neelamegam & Shankar (2015) (blue) and Chandra et al. (2019) (red). The different black symbols represent the transition Reynolds number from friction factor $(\circ )$, dye-stream visualisation $(\triangle )$ and detection of wall oscillations $(\boldsymbol {\nabla })$. Redrawn and augmented, with permission, from Verma & Kumaran (2012).

A global stability analysis for the flow in a deformed tube has not been carried out so far. An approximate analysis for a quasi-fully developed flow has been carried out by Verma & Kumaran (2015). This is valid when the slope of the wall $\alpha$ is small, so that the flow can be considered locally parallel. Even when the slope is small, the velocity profile could be very different from the parabolic profile for ${Re} \alpha \sim 1$, because the inertial terms in the momentum equations scale as ${Re} \alpha$. The detailed explanation of this procedure is provided for a channel flow in the next subsection, and so it is not repeated here.

The stability analysis of Verma & Kumaran (2015) predicted that the growth rate of perturbations first becomes positive in the downstream converging end of the test section; this is in agreement with the visual observations in, for instance, figure 12. The value of the transition Reynolds number was in quantitative agreement with the experimental results. It was also reported that there is agreement only if the effect of tube deformation on the velocity profile and pressure gradient is considered in the linear stability analysis; for a parabolic flow and a constant pressure gradient, the transition Reynolds number is higher than the experimental value by a factor of 10. An analysis of the eigenfunctions of the most unstable mode in Verma & Kumaran (2015) indicated that the perturbations are confined to thin regions both in the fluid and in the gel, suggesting that transition is due to the wall-mode instability.

8.2.2. Channel flow

The transition in a microchannel of height (smallest dimension) $100$ and $160\ \mathrm {\mu }$m, width about 1.5 mm and length about 3 cm with three hard walls and one soft wall was studied by Verma & Kumaran (2013). An indentation of depth 100 or 160 $\mathrm {\mu }$m in the shape of a channel was made in a hard PDMS stamp of shear modulus 0.55 MPa; one of the channel shapes, the two-inlet shape, is shown in figure 15(a). The stamp made of hard PDMS was bonded with a layer of soft PDMS of the desired shear modulus mounted on a glass slide. This results in a rectangular channel of height $100 - 160\ \mathrm {\mu }$m and width 1.5 mm in a gel block, three walls of which are made of hard gel and one of soft gel. A schematic of the cross-section is shown in figure 15(b). The top view of the set-up is shown in figure 15(c), where the microchannel within the gel block is highlighted by blue ink.

Figure 15. (a) Top view of the microchannel configuration with Y inlet. (b) Cross-section showing three hard walls and one soft wall. (c) Image from the top of the gel block with blue dye showing the location of the microchannel. Redrawn with permission from Verma & Kumaran (2013).

The microchannel has a development section of length about 0.8 cm, where all four walls are made of hard gel, to ensure that disturbances are damped out before the fluid enters the test section. The development section is fabricated by adding extra catalyst to a portion of the soft PDMS film before it is bonded to the hard PDMS stamp. The microchannel also contains a pressure port at the entrance to the test section, which can be used for the inlet pressure measurement. The channel height, which is $h_0 = 100/160\ \mathrm {\mu }$m in the absence of flow, can increase in the test section when a pressure difference is applied across the ends. However, the Reynolds number based on the flow rate and the channel width $W$, ${Re} = (\rho _f Q / \mu _f W)$, does not change when there is a change in the channel height.

Transition was detected in two ways. The first used a dye stream, where black dye is injected through one inlet in figure 15(a) and clear water through the other inlet. Images taken from above are used to examine whether there is cross-stream mixing across the microchannel. The second method is the measurement of the pressure drop as a function of the flow rate. When there is substantial wall deformation, friction factor relations derived for a rigid channel cannot be used to indicate transition; as discussed below, a more sophisticated procedure is required. In order to detect wall oscillations, fluorescent microbeads are embedded in the soft wall of the microchannel. These are illuminated from the side by a laser beam, and light scattered from the microbeads is imaged from above. The onset of wall oscillations is indicated by a sharp increase in the root mean square of the intensity fluctuations of the scattered light.

Figure 16 shows images, from above, of the flow in the microchannel at different flow rates when the soft wall is made of PDMS gel with shear modulus 18 kPa. Water with black dye and clear water are pumped into the lower and upper inlets at equal flow rates. The images are centred at the locations L1, L2, L3 and L4 in figure 15(a). Note that the location L1 is in the development section, where all four walls are made of hard gel, while the other locations are in the test section. The edges of the microchannel are shown by the black dashed lines, and the width of 1.5 mm provides a scale for the images. At a low flow rate, it is observed that there is no cross-stream mixing over the time required for the fluid to flow through the microchannel. As the flow rate is increased, there is the onset of intense cross-stream mixing at a Reynolds number of about 200 at the downstream locations in the channel. Even when there is near-complete mixing downstream, there is virtually no mixing at the upstream location L2 at ${Re} = 200$. When the Reynolds number is increased, intense mixing is observed at the upstream location L2 as well, but there is no mixing in the development section, indicating that the flow is laminar. This intense cross-stream mixing at a low Reynolds number of 200 indicates a transition from a laminar flow to a more complicated flow profile; the transition Reynolds number is five times lower than the value of 1000 for a rigid channel. The transition Reynolds number was found to depend on the shear modulus of the soft wall – it was about 289 when the wall had a shear modulus of about 28 kPa and 311 when the wall had a shear modulus of about 35 kPa.

Figure 16. Top view of the mixing between a clear water stream and a stream of water coloured with black dye in the Y-channel configuration shown in figure 15(a) when the compliant wall is made of gel with shear modulus about 18 kPa. Reprinted, with permission, from Verma & Kumaran (2013).

There is substantial deformation of the soft wall due to the applied pressure gradient, as shown in figure 17. Figure 17(a) shows the side view of the channel in the absence of flow, while figure 17(b) shows the side view when the Reynolds number is 222 in a channel in which the soft wall is made of material of shear modulus about 18 kPa. It is observed that there is virtually no expansion in the development section, but there is significant expansion in the test section. The channel first diverges at the entrance to the test section, and then converges further downstream as the fluid pressure decreases. The maximum height of the channel is about 3–4 times the height of the undeformed channel in this particular case. However, the expansion is only in the soft wall, and the other three walls are undeformed.

Figure 17. Side view of the channel with soft wall made of material of shear modulus 18 kPa in the absence of flow (a) and the deformed channel in the presence of flow at a Reynolds number of 200 (b); the side view (c) and cross-section at different downstream locations (d) of the deformed channel reconstructed in a simulation. Reprinted, with permission, from Verma & Kumaran (2013).

Due to the substantial deformation, the laminar flow in the channel is not a Poiseuille flow, and it is necessary to compute the velocity profile and the pressure drop for the deformed channel. The channel is reconstructed from the dimensions measured in the side view in figure 17(b), and the reconstruction assumes that the three rigid walls are undeformed and only the soft wall is deformed. The side view of the reconstructed channel is shown in figure 17(c) and the cross-sections at different downstream locations are shown in figure 17(d). The reconstructed shape of the channel was used in computational fluid dynamics (CFD) simulations employing the ANSYS FLUENT software in order to determine the velocity profile and the pressure as a function of the downstream location. In the CFD simulations, the flow is assumed to be laminar, a constant flow rate boundary condition with a plug flow is applied at the inlet and zero pressure is applied at the outlet. The velocity profiles are different from a parabolic profile in a channel with parallel walls – they have a gentler gradient at the wall and a higher curvature at the centre in the upstream diverging section, and a sharper gradient at the wall and a lower curvature at the centre in the downstream converging section. The pressure is not a linear function of downstream distance; the pressure gradient is small in the upstream expanding section of the channel and it increases in the downstream converging section.

Since the pressure difference across the ends of the channel is calculated assuming the flow is laminar in the CFD simulations, the experimental pressure drop will be higher than that predicted for a laminar flow after transition. Therefore, the transition Reynolds number based on the pressure difference is the Reynolds number at which the experimental pressure drop is higher than that predicted by the laminar flow simulations for a channel of the same shape. In addition to the dye-mixing method and the detection of wall oscillations, this provides the third method for detecting transition used by Verma & Kumaran (2013). Verma & Kumaran (2013) also carried out a linear stability analysis using the parallel flow approximation locally at different streamwise locations along the channel. The linear analysis assumes that the flow is two-dimensional and fully developed with a mean velocity profile and pressure gradient obtained, at each location, from the CFD simulations for a laminar flow. This is valid only if the wavelength of the perturbations is much smaller than the length scale for flow development and the spanwise extent of the channel. In figure 18(a), the experimental results for the transition Reynolds number from the dye-mixing method and the pressure difference method are compared with the theoretical prediction from the linear stability analysis. The results from the three methods are in quantitative agreement. From the linear stability analysis, it was found that the flow first goes unstable at the downstream converging section of the microchannel. The Reynolds numbers for the onset of an instability at the outlet ($x = 3$ cm) and at a location 3 mm upstream of the outlet ($x = 2.7$ cm) are superposed on the experimental results. The experimental results for the transition Reynolds number are bracketed by results of the linear stability analysis at the two downstream locations. This indicates that the linear stability analysis has some predictive capability for predicting the transition in the downstream converging section.

Figure 18. The transition Reynolds number as a function of the parameter $\varSigma$ for the flow in a microchannel. (a) The transition Reynolds number measured by different experimental methods superposed on the theoretical predictions at the downstream locations $x = 3$ cm and $x = 2.7$ cm shown in figure 17(c) from Verma & Kumaran (2013). (b) The transition Reynolds number from Verma & Kumaran (2013) (blue), Srinivas & Kumaran (2017b) (black) and Kumaran & Bandaru (2016) (red). The horizontal dashed line is the hard-wall laminar–turbulent transition Reynolds number and the inclined dashed line is the scaling law ${Re}_t \propto \varSigma ^{5/8}$. Panel (a) redrawn and augmented with permission from Verma & Kumaran (2013).

The transition Reynolds numbers from the study of Verma & Kumaran (2013) are plotted along with those of Srinivas & Kumaran (2017b) for larger channels of height 0.6–1.8 mm made of polyacrylamide, and smaller channels of height 35–100 $\mathrm {\mu }$m in the study of Kumaran & Bandaru (2016), in figure 18(b). These studies cover a range of about three orders of magnitude in the parameter $\varSigma$ and about two orders of magnitude in the transition Reynolds number. The scaling law ${Re}_t \propto \varSigma ^{5/8}$ appears to fit the data reasonably well over this entire range, though there is a spread in the data. One reason for the spread is that the experiments are carried out at discrete values of the Reynolds number separated by 50–100, and so the precise transition Reynolds number is difficult to pinpoint. The other is the variation in the deformed channel configuration in different experiments at the same Reynolds number – the deformed configuration depends on the cross-section, the length and height of the channel and the elasticity of the material used.

An interesting observation in figure 18(b) is the presence of transition even above the hard-wall laminar–turbulent transition Reynolds number of 1000 in the experiments of Srinivas & Kumaran (2017b) for the two experiments with the highest $\varSigma$. It is observed that there is a transition between two distinct types of turbulence, hard-wall and soft-wall turbulence, at a higher Reynolds number compared with the hard-wall laminar–turbulent transition. This is discussed in § 8.3.

8.3. Turbulence

There have been only two studies, those of Srinivas & Kumaran (2015) and Srinivas & Kumaran (2017b), on the nature of the flow after transition. Both of these studies indicate that the flow can be characterised as a turbulent flow, but there are qualitative differences from the turbulent flow in a rigid channel.

The analysis of a turbulent flow by Srinivas & Kumaran (2015) was carried out for the same set-up as Verma & Kumaran (2013). The fluid velocities were measured using particle image velocimetry (PIV) with polymer microspheres of diameter about 2 $\mathrm {\mu }$m as seed particles for flow visualisation. Recall that the channel has one soft wall at the bottom and three hard walls, and has a rectangular cross-section with width 1.5 mm and height 160 $\mathrm {\mu }$m in the absence of flow. A laser sheet from above was used to illuminate the channel, as shown in figure 19(b), and masking tapes were used to restrict the laser illumination to a width of 200 $\mathrm {\mu }$m along the mid-plane of the channel in the spanwise direction. A nanosecond pulsed laser was used to generate pairs of images separated by a time interval of a few nanoseconds. Two-dimensional images were taken from the side through a high-resolution camera at three different locations in the streamwise direction shown in figure 19(a), and these were used to generate the velocity profiles. In the velocity profiles, the height across the microchannel is measured from the bottom (soft) surface, as shown in figure 19(b).

Figure 19. The top view (a) and side view (b) of the microchannel and the PIV camera and laser sheet for the measurement of turbulence statistics in a microchannel. The dimensions of the cross-section are the same as those in figure 15(b). Redrawn with permission from Srinivas & Kumaran (2015).

The salient features of the turbulent flow are summarised in figure 20, which shows the turbulent statistics for a flow through a channel with one soft wall made of PDMS gel with elasticity modulus 18 kPa at the downstream location C in figure 19(a). These results, redrawn from Srinivas & Kumaran (2015), are shown for three Reynolds numbers based on the flow rate and channel width, one (${Re} = 222$) lower and the other two (${Re} = 277$ and $415$) higher than the transition Reynolds number. The location $y=0$ is the soft wall in figure 19(b), and the location of the top hard wall is shown by the dashed lines on the right. The latter changes due to increased channel deformation as the Reynolds number is increased. In all cases, the symbols are the experimental results, and the error bars show one standard deviation over three independent runs carried out on the same channel geometry for the same flow rate.

Figure 20. The mean velocity (a), mean square of the streamwise fluctuating velocity (b), the Reynolds stress (c) and the turbulent energy production rate (d) obtained from PIV measurements along the centreline of the microchannel at location C in figure 19(b), for a laminar flow at Reynolds number $222$ ($\circ$, blue) and for turbulent flows at Reynolds number $277$ ($\triangle$, red) and $415$ ($\boldsymbol {\nabla }$). Redrawn with permission from Srinivas & Kumaran (2015).

The mean velocity profiles are shown by the symbols in figure 20(a). The solid lines are the velocity profiles obtained from simulations for the laminar flow at the same location for the flow in a channel of the same shape as that in the experiments. These were obtained as explained in the discussion around figure 17. It is observed that there is quantitative agreement between the experimental and simulation mean velocity profiles for ${Re} = 222$. However, after transition, the experimental result is not in agreement with the laminar velocity profile. The experimental result is more plug-like, with a smaller curvature in the centre and larger gradients near the walls in comparison to the laminar flow – the experimental profile is similar to the velocity profile for a turbulent flow in a rigid channel.

Figures 20(b) shows the root mean square of the velocity fluctuations in the streamwise direction. It should be noted that the intensity of the fluctuations is measurable even for the pre-transition flow at ${Re}=222$ – as explained at the beginning of this section, flow conditioning is not perfect in a channel of submillimetre size with compliant walls. However, there is a clear increase in the intensity of the fluctuations after transition at ${Re} = 277$, and a further increase at ${Re} = 415$. The near-wall maximum in the streamwise root mean square velocity, which is a characteristic of turbulent channel flows, is clearly observed in figure 20(b). More interestingly, the profiles of the velocity fluctuations are not symmetric, and the maximum near the bottom (soft) wall is about two times larger than that near the top (hard) wall. This asymmetry is observed in the profiles of the cross-stream root mean square velocity as well. This clearly indicates that the dynamics of the bottom wall plays a significant role in the generation of turbulent fluctuations. In figure 20(b), the velocity statistics could not be obtained within a region of thickness about $10\ \mathrm {\mu }$m from the wall, due to limitations with the resolution of the measurements and the size of the seed particles. Therefore, it is not clear if the root mean square velocity decreases to zero at the wall.

The correlation $\overline {v_x' v_y'}$, which is the ratio of the Reynolds stress and the fluid density, is shown as a function of the cross-stream distance in figure 20(c). The Reynolds stress is zero, to within the experimental error bars, for the laminar flow at ${Re} = 222$. After transition, the characteristic near-wall maximum of the Reynolds stress is observed. The Reynolds stress profiles are asymmetric, with the stress near the bottom (soft) wall significantly higher than that near the top (hard) wall. The Reynolds stress does decrease to zero at the top wall. However, it is difficult to extrapolate the Reynolds stress to zero at the bottom (soft) wall, though it is not possible to rule out a sharp decrease in the Reynolds stress to zero. Figure 20(d) shows the turbulent energy production rate, which is the product of the Reynolds stress and the strain rate. If the Reynolds stress is non-zero at the bottom wall, the turbulent energy production rate is a maximum at the bottom wall. This is in contrast to the flow past a rigid surface, where the turbulent energy production rate has a near-wall maximum and it decreases to zero at the wall.

The turbulent flow past a rigid surface contains a logarithmic layer for $5 < y_+ < 30$, where the velocity profiles are described by the logarithmic law ($\bar {v}_x/v \approx 2.5 \log {(y v_\ast / \nu )} + 5$). This law is found to apply at Reynolds number greater than about 3500. For the flow in a soft-walled channel, Srinivas & Kumaran (2015) reported that the flow does exhibit a logarithmic layer close to the wall, but the logarithmic law is different from that in a hard channel, and it does depend on the wall shear modulus. This is an aspect that is intriguing, but which requires a lot more confirmation.

In summary, the flow after transition in a soft-walled microchannel at Reynolds number as low as 250 does exhibit many of the features of a turbulent flow in a rigid channel, including the plug-like nature of the mean velocity profile, the near-wall maximum in the streamwise root mean square velocity and the characteristic shape of the Reynolds stress curves. However, it is clear that wall flexibility does play a role in generating turbulence, because the velocity fluctuations are significantly higher near the soft wall in comparison with the rigid surface. There are still unanswered questions, including whether the magnitude of the root mean square velocity does decrease to zero at the flexible surface, and whether there is a viscous sublayer at the wall. However, it is clear that that the flow after transition can be accurately characterised as turbulent, but one that is qualitatively different from the turbulent flow in a rigid channel.

The magnitudes of the velocity fluctuations in this turbulent flow, when scaled by the average velocity, are rather large. The background level of fluctuations in a soft-walled channel is rather large even for a laminar flow. After transition, the level of the fluid velocity fluctuations is much larger than that in a rigid channel. At Reynolds number 250–400 in a microchannel, when scaled by the suitable powers of the average velocity, the mean square velocities are comparable to or larger than those in a rigid channel at Reynolds number in the range 12 000–20 000.

Srinivas & Kumaran (2017b) carried out experiments in a larger rectangular channel of height 0.6 and 1.8 mm and width 1.3 cm made of polyacrylamide. The cross-section and the dimensions of the channel are shown in figure 21. There was a development section of length about 13 cm where the walls were made of hard gel of shear modulus about 16 kPa, followed by a test section where the gel was made with shear modulus 0.75 or 2.18 kPa. The velocity was measured using PIV with a set-up similar to that in figure 19, using glass beads of diameter 8–14 $\mathrm {\mu }$m as seed particles. Compared with the flow in a microchannel, a much larger Reynolds number could be examined in the channel of larger dimension, and the wall displacement was much larger due to the lower shear modulus of polyacrylamide. The range of Reynolds numbers in the experiments exceeded the Reynolds number of 1000 for the laminar–turbulent transition in a rigid channel, and so the effect of wall compliance on the transition in a rigid channel could also be examined. In contrast to the microchannel geometry in figure 19, all four walls of the channel are made with the same shear modulus. However, there was an asymmetry between the top and bottom walls because the bottom wall was mounted on a rigid glass plate, while the top wall was unrestrained.

Figure 21. Channel top view (a), cross-section (b) and deformed cross-section (c) used in the experiments of Srinivas & Kumaran (2017b). Redrawn with permission.

The average and the root mean square of the fluctuations in the wall displacement were measured by optical imaging in the direction parallel to the surface using a camera looking vertically downward. Glass beads of diameter 8–14 $\mathrm {\mu }$m were embedded in the gel wall, and the tangential displacement was measured separately on the top and bottom surfaces by focusing the camera appropriately. The displacement perpendicular to the surface was measured on the top and bottom surfaces using a camera from the side. The results for the root mean square of the displacement fluctuations are shown in figure 22 for a channel with undeformed height 1.8 mm and made of material of shear modulus 0.75 kPa in the test section. The tangential displacement is shown in figure 22(a) separately for the top and bottom walls. The transition Reynolds number ${\sim }1000$ for a rigid tube is shown by the line labelled HW in figure 22. There is no evident increase in the displacement fluctuations in either the tangential or the normal direction at this Reynolds number, although, as discussed next, the fluid velocity fluctuations do exhibit a discontinuous change. Figure 22(a) shows that there is an increase in the streamwise displacement fluctuations at a higher Reynolds number of about 1400, labelled SW. The increase is approximately the same on the top and bottom walls. However, there is no discernible cross-stream displacement fluctuation in figure 22(b). When the Reynolds number is increased beyond 1700, labelled WF (wall flutter), there is an increase in the displacement fluctuations both parallel and perpendicular to the surface at the top wall. There is no commensurate increase in the displacement fluctuation at the bottom wall. Srinivas & Kumaran (2017b) reported that there are visible travelling waves on the top wall alone. The asymmetry between the top and bottom wall is because the bottom wall is fixed to a glass substrate, while the top wall is unrestrained; when the top wall is also constrained with a glass plate, the wall flutter is not observed.

Figure 22. The root mean square of the displacement parallel (a) and perpendicular (b) to the surface in a channel of height 1.8 mm in the absence of flow at downstream location III in figure 21(a) when the wall is made of polyacrylamide gel of shear modulus 0.75 kPa. In (a), the solid and dashed lines are the tangential displacements on the top and bottom walls, and (b) shows the normal displacement on the top wall. Redrawn with permission from Srinivas & Kumaran (2017b).

The velocity statistics also show signs of two different transitions, after the laminar–turbulent transition at a Reynolds number of about 1000, as shown in figure 23. The flow is laminar at a Reynolds number of 768 and there is a transition to a turbulent flow when the Reynolds number is increased to 1071. Here, the mean velocity is distinctly different from the parabolic profile, shown by the dashed line in figure 23(a). There is a significant increase in the velocity fluctuations in the streamwise and cross-stream directions, and the large near-wall maximum in the streamwise velocity fluctuations is clearly visible in figure 23(b). At this Reynolds number, there are no discernible wall fluctuations in figure 22. When the Reynolds number is increased from 1332 to 1515, the onset of tangential wall oscillations in figure 22 is accompanied by a change in the shape of the mean velocity profile in figure 23(a), a significant increase in the streamwise velocity fluctuations in figure 23(b) and a significant increase in the Reynolds stress in figure 23(d). However, the profiles are still symmetric about the centreline of the channel. When the Reynolds number is further increased from 1734 to 1973, there is a significant expansion of the channel accompanying the onset of normal wall oscillations in figure 22(b). There is a visible asymmetry in the profiles for the mean and fluctuating velocities and in the Reynolds stress. The velocity fluctuations near the top wall are significantly higher than those near the bottom wall; this is due to the presence of significantly enhanced wall displacement fluctuations at the top wall.

Figure 23. The mean velocity (a), root mean square of the streamwise (b) and cross-stream (c) fluctuating velocities and the Reynolds stress (d) at Reynolds number $Re=768$ ($\circ$), $Re=1071$ ($\triangle$, brown), $Re=1332$ ($\boldsymbol {\nabla }$, brown), $Re=1515$ ($\triangleleft$, blue), $Re=1734$ ($\triangleright$, blue) and $Re=1973$ ($\diamond$, red) in a channel of height 1.8 mm in the absence of flow at downstream location III in figure 21(a) when the wall is made of polyacrylamide gel of shear modulus 0.75 kPa. Redrawn with permission from Srinivas & Kumaran (2017b).

The study of Srinivas & Kumaran (2017b) shows that there are at least two distinct transitions in the turbulent flow in a channel in addition to the hard-wall laminar–turbulent transition at ${Re} = 1000$. The transitions occurred at a Reynolds number lower than $1000$ for the channel with height $0.6$ mm and at a Reynolds number higher than $1000$ for the channel with height $1.8$ mm. The first transition, called the ‘soft-wall’ transition, has the characteristics of the wall mode instability, where there are visible displacement fluctuations in the tangential direction, but no discernible fluctuations perpendicular to the surface. This is accompanied by a sharp increase in the velocity fluctuations in the fluid and in the Reynolds stress, but the profiles are symmetric about the centreline. There is onset of flutter on the top surface at the second ‘wall flutter’ transition, and this occurs only if the top surface is unrestrained. The significant observation in Srinivas & Kumaran (2017b) was that the soft-wall transition could take place even in a turbulent flow after the hard-wall laminar–turbulent transition. The coupling between the wall and fluid dynamics after transition qualitatively alters the nature of the turbulent flow, and the post-transition flow, which exhibits features distinct from the turbulence in the flow past rigid surfaces, was characterised as ‘soft-wall turbulence’.

8.4. Viscoelastic fluids

The flow of viscoelastic fluid past a polymer gel in the viscous limit was studied by Neelamegam, Shankar & Das (2013) in the context of the ‘purely elastic’ instability due to the hoop stress in a rotating polymeric fluid (Shaqfeh 1996; Groisman & Steinberg 2004). When an elastic fluid is placed between two rigid plates as shown in figure 9(a) and the strain rate is increased, the flow transitions from a simple viscometric flow through a series of bifurcations to a state called ‘elastic turbulence’. This is indicated by sharp increases in the apparent viscosity for specific values of the Weissenberg number, which is the product of the strain rate and the polymer relaxation time. When the bottom rigid plate is replaced by a soft elastic gel, as shown in figure 9(b), Neelamegam et al. (2013) found that the elastic instability is suppressed, and the sharp increase in the apparent viscosity is absent. At fixed strain rate, the presence of a soft bottom surface also suppresses temporal fluctuations in the stress which are observed when a rigid bottom surface is used. Using flow visualisation, secondary flow patterns observed in the flow between two rigid surfaces were absent when one of the surfaces is a soft gel. An important difference between a Newtonian and a polymeric fluid is the first normal stress difference, which results in an outward force on the top and bottom surface. The mechanism for the suppression of the instability suggested by the authors is as follows. When the bottom surface is made of soft gel, the normal stress difference could compress the gel, increase the width of the fluid layer and thereby decrease the actual shear rate. Due to this, the actual Weissenberg number is lower than that calculated assuming that the width of the fluid layer is unchanged, and therefore the instability is not observed.

The effect of polymers on the transition and turbulence in a compliant channel was studied by Srinivas & Kumaran (2017a) for the flow through a PDMS microchannel shown in figure 15 and by Chandra et al. (2019) for the flow in a PDMS microtube of diameter 300–400 $\mathrm {\mu }$m. The elasticity modulus of the PDMS gel was 0.55 MPa for the hard walls and about 18 kPa for the soft walls in Srinivas & Kumaran (2017a), and in the range 30–50 kPa in Chandra et al. (2019). The viscoelastic fluid consisted of polyacrylamide polymers (not cross-linked) of molecular weight $4 \times 10^{4}$ and $5 \times 10^{6}$ dissolved in water. The concentration of the polymer was very small. In Srinivas & Kumaran (2017a), concentrations up to 1500 and 50 p.p.m. were used for the polymer with molecular weight $4 \times 10^{4}$ and $5 \times 10^{6}$, respectively, and the steady shear viscosity did not exceed the viscosity of water by more than 10 %. Chandra et al. (2019) used concentrations up to 2500 p.p.m. for polyacrylamide with molecular weight $5 \times 10^{6}$, and the maximum viscosity was about five times the viscosity of water.

An important material parameter is the polymer relaxation time. There are different methods to measure the polymer relaxation time, such as small-amplitude oscillatory measurements fitted to the Maxwell model (Doi & Edwards 1988), first normal stress difference (Larson, Shaqfeh & Muller 1990) and capillary break-up (CaBeR) technique for the extensional relaxation time (Dinic et al. 2015). Zell et al. (2010) found that there are significant differences in the results of different measurements. Srinivas & Kumaran (2017a) used the Zimm model based on the polymer radius of gyration and the solvent viscosity in order to estimate the relaxation time and found that the relaxation time is about 30 ms for polyacrylamide with molecular weight $5 \times 10^{6}$ and $5\ \mathrm {\mu }$s for polyacrylamide with molecular weight $4 \times 10^{4}$. Chandra et al. (2019) used small-amplitude oscillatory measurements fitted to the Maxwell model to determine the relaxation time and found a concentration-dependent relaxation time ranging from 1.4 ms for a 100 p.p.m. solution to about 7 ms for a 2500 p.p.m. solution for polyacrylamide with molecular weight $5 \times 10^{6}$. The variation in relaxation time also results in a variation in the the elasticity number, ${El} = (\lambda \mu _f /\rho _f h_f^{2})$, where $\lambda$ is the polymer relaxation time. The elasticity number in the experiments of Srinivas & Kumaran (2017a) was in the range 0.12–1.2 for polyacrylamide with molecular weight $5 \times 10^{6}$ and $2 \times 10^{-5} - 2 \times 10^{-4}$ for polyacrylamide with molecular weight $4 \times 10^{4}$. In the experiments of Chandra et al. (2019), the elasticity number varied in the range $8 \times 10^{-3}$ for a 50 p.p.m. solution to 0.25 for a 2500 p.p.m. solution. Thus, the elasticity number has been varied over a very wide range in the experiments conducted so far. Since the solution viscosity considered by Chandra et al. (2019) varied over a wider range from $1\times 10^{-3}$ to $7 \times 10^{-3}$ kg m$^{-1}$ s$^{-1}$, the parameter $\beta$ in the Oldroyd model, which is the ratio of the solvent viscosity and the solution viscosity, was also varied over a wide range. Srinivas & Kumaran (2017a) considered a very dilute solution, so $\beta$ was close to $1$.

An important parameter is the ratio of the polymer concentration and the ‘overlap’ concentration $c^{\ast }$ (de Gennes 1976, 1979). In a dilute solution with concentration below the overlap concentration, a polymer molecule in the solution does not interact with other molecules, whereas interactions between polymer molecules are important in a semi-dilute solution when the concentration exceeds the overlap concentration. The overlap concentration was estimated as 150 p.p.m. for polyacrylamide with molecular weight $5 \times 10^{6}$ and about 1670 p.p.m. for polyacrylamide with molecular weight $4 \times 10^{4}$. Therefore the experiments of Srinivas & Kumaran (2017a) were entirely in the dilute regime, whereas the experiments of Chandra et al. (2019) extended over the dilute and semi-dilute regimes.

The experimental set-up of Srinivas & Kumaran (2017a) was the same as that of Verma & Kumaran (2013) shown in figure 15. Transition was detected by dye-mixing experiments of the type shown in figure 16, as well as the velocity statistics measured using PIV as shown in figure 19. Srinivas & Kumaran (2017a) found that the addition of polymer has no effect on the transition Reynolds number if the concentration is below a threshold value, which is about 1 and 500 p.p.m. for polyacrylamide with molecular weight $5 \times 10^{6}$ and $4 \times 10^{4}$, respectively. Above this threshold concentration, there is a systematic decrease in the transition Reynolds number. For the highest concentration of 1500 p.p.m. for the polymer with molecular weight $4 \times 10^{4}$, the transition Reynolds number is 139, while for the highest concentration of 50 p.p.m. for the polymer with molecular weight $5 \times 10^{6}$, the transition Reynolds number is 39. This is in contrast to the transition Reynolds number of 291 for the same channel for a Newtonian fluid. Thus, addition of small amounts of polymer results in a decrease of an order of magnitude in the transition Reynolds number. This was consistent with the prediction of Chokshi et al. (2015) for the effect of viscoelasticity on the wall mode instability in a channel.

The second important observation in Srinivas & Kumaran (2017a) was that the polymer has a significant damping effect on the fluid velocity fluctuations. Figure 20 shows that there are significant fluctuations even in a laminar flow. Even when the polymer concentration is below the threshold concentration and there is no reduction in the transition Reynolds number, there is a reduction by a factor of $2$ in the root mean square of the velocity fluctuations in comparison to the laminar flow of pure water. In the turbulent regime, there is a reduction by a factor of 5 in the root mean square velocities and by a factor of 10 in the Reynolds stress. The equivalent of the images shown in figure 15 for the flow of polymer solutions do not exhibit the intense cross-stream mixing across the span, and the instability is indicated by waves at the interface between the two streams with and without dye. Chandra et al. (2019) detected transition in the tube flow based on the departure of the velocity profile from the parabolic profile, as well as the intensity of the velocity fluctuations scaled by the average velocity. For dilute polymer solutions with concentration less than about 100 p.p.m., there was a reduction in the transition Reynolds number due to the addition of the polymer solution. The scaling of the transition Reynolds number was ${Re}_t \propto \varSigma ^{0.7}$ for a viscoelastic fluid, which is not very different from the scaling ${Re}_t \propto \varSigma ^{0.65}$ for a Newtonian fluid and the scaling ${Re}_t \propto \varSigma ^{5/8}$ found by Verma & Kumaran (2012). Chandra et al. (2019) concluded that the destabilisation is due the wall mode instability modified by fluid viscoelasticity, as reported by Srinivas & Kumaran (2017a). For higher polymer concentration, Chandra et al. (2019) found a very different cause of destabilisation, which is the elasto-inertial instability (Zakin et al. 1977; Samanta et al. 2013) due to combined effects of elasticity and inertia. Due to this instability, the transition Reynolds number for the flow of a viscoelastic fluid through a rigid tube could be lower than that for a Newtonian fluid. The transition Reynolds number follows the scaling ${Re}_t \propto {El} (1-\beta )^{-1/2}$ (Chandra, Shankar & Das 2018), where ${El}$ is the elasticity number and $\beta$ is the ratio of the solvent and solution viscosity. For the flow in a compliant tube, when the polymer concentration exceeded 150 p.p.m., Chandra et al. (2019) found the same scaling as that for a rigid tube, but the transition Reynolds number was lower than that for a rigid tube. This indicates that the wall flexibility destabilises the elasto-inertial mode, resulting in a lower transition Reynolds number.

8.5. Mixing

Shrivastava, Cussler & Kumar (2008) carried out an interesting study of the mass transfer enhancement in the pressure-driven flow in a channel in which one wall was made of soft PDMS gel with shear modulus in the range 1–3 kPa. The channel of length 10 cm, width 4 cm and height 1.5 mm was made using polycarbonate blocks, and the working and counter electrodes of dimension $3\ \textrm {cm} \times 1\ \textrm {cm}$ made of platinum foil of thickness 0.05 mm were glued to the top and bottom walls. The gel was placed on the counter electrode, and a glycerol-based electrolyte solution was pumped through the channel. The Reynolds number was in the range 0.1–60, and the dimensionless parameter $\varGamma = (\mu _f V_f/G h_f)$ was in the range $10^{-3} - 10$. A voltage sweep was applied across the reference and working electrodes to determine the limiting current, and this was related to the mass transfer coefficient at the gel surface. Control experiments were also carried out where the gel block was replaced by a rigid block. For $\varGamma < 0.1$, the authors found that the limiting current is the same for the flow past gels and rigid blocks, but as $\varGamma$ increases, the mass transfer coefficient past a gel surface is significantly enhanced in comparison to that past a rigid block. The ratio of the currents in the flow past gel and rigid blocks was plotted as a function of $\varGamma$. It was found that the ratio is $1$ for small values of $\varGamma$, but it increases continuously to about $1.25$ for the highest value of $\varGamma = 10$. This indicates that there is enhanced mass transfer in the flow past a gel surface in comparison to the flow past a rigid surface. Though the authors attributed this to an instability at the fluid–gel interface, they did not report any discontinuous change in the mass transfer coefficient at a critical value of $\varGamma$. Nevertheless, this was the first study to suggest the possibility of using a soft gel to enhance mass transfer from a surface in small-scale flows.

Studies of cross-stream mixing across a microchannel have been carried out in Y-shaped channels of the kind shown in figure 15. Image analysis has been used on images of dye-mixing experiments (figure 16, for example) to quantify mixing. Images at different locations, shown by the red rectangles in the top panel in figure 16, are captured and the mean and standard deviation of the pixel intensity are calculated. The extent of mixing is related to the standard deviation of the pixel intensity about the mean – a lower standard deviation indicates less variation in the pixel intensity across the channel, and therefore better mixing. Verma & Kumaran (2013) defined a segregation index based on image analysis, and showed that the segregation index decreases discontinuously at transition.

Another method is to pump two contrasting fluids through the two inlets in figure 15, and to analyse the fluids emerging from the two outlets separately. The fluids at the two outlets will be a homogeneous mixture of the inlet fluids if the mixing is perfect, whereas the fluid at the outlet will be the same as the fluid at the corresponding inlet when there is no mixing. Verma & Kumaran (2012) pumped in deionised water in one inlet and a solution of tannic acid in the other inlet, and determined the concentration of tannic acid at the outlet from conductivity measurements. Kumaran & Bandaru (2016) carried out experiments in smaller microchannels of length 2 cm, width 0.5 mm and height 35–100 $\mathrm {\mu }$m with two inlets and two symmetric outlets, in which one soft wall was made with shear modulus 17 kPa. The two contrasting fluids used in Kumaran & Bandaru (2016) included a combination of one acid and one base where the pH at the outlet was measured, and a combination of a suspension of $6\ \mathrm {\mu }$m polystyrene particles and clear fluid where the concentration of the particles was measured at the outlet. The authors reported that for the smallest microchannel of height 35 $\mathrm {\mu }$m, there is no cross-stream mixing when the Reynolds number is 66, but there is imperfect cross-stream mixing when the Reynolds number is 133, and there is perfect cross-stream mixing (the two fluids are identical to within experimental resolution) at a Reynolds number as low as 200.

A striking conclusion is that this ultrafast mixing across the 0.5–1.5 mm width of the microchannel takes place within a time interval of a few tens of milliseconds. The flow velocity in the microchannels is of the order of 1 m s$^{-1}$, while the length of the channel is 2–3 cm. Therefore, the residence time of the fluid in the channel is of the order of 20–30 ms. Within this time period, there is near-perfect mixing across the width of the microchannel. This is in contrast to mixing by molecular diffusion in a laminar flow which is a very slow process. A characteristic molecular diffusion coefficient for solutes with small molecular weight in water is $10^{-9}$ m$^{2}$ s$^{-1}$. The characteristic time for molecular diffusion across a width of 1 mm is the ratio of the square of the width and the diffusion coefficient, which is $10^{3}$ s. Thus, the time for mixing due to the soft-wall turbulence is five orders of magnitude less than that for diffusion in a laminar flow.

This enhanced mixing requires a less-than-proportionate increase in the pumping cost due to the high flow rate. For a rigid channel, the pressure difference required for driving the flow is proportional to the flow rate. When the channel is made of a compliant surface, there is expansion of the channel, as shown in figure 17. Due to this, there is an increase in the average area of cross-section and a decrease in the pressure difference. Thus, using a soft wall to induce turbulence is a simple and low-cost method for achieving cross-stream mixing in microfluidic devices, where slow mixing due to molecular diffusion is a technological barrier for sample preparation. In Abbasi et al. (2019), the utility of this technology has recently been demonstrated for a point-of-care device.

9. Simulation of turbulent flows

The first direct numerical simulation of the flow in a channel with a compliant wall was carried out by Xu, Rempfer & Lumley (2003) and Rempfer et al. (2003). The authors used a pseudospectral method for solving the equations that is commonly used for direct numerical simulations (Moin, Kim & Moser 1987) involving Fourier transforms in the periodic streamwise and spanwise directions and a Chebyshev transform in the wall-normal direction. These transforms can be carried out only if the domain is cuboidal, and it cannot be applied in an irregular domain formed by a deformed wall. Therefore, the authors used a mapping to convert the irregular domain into a cuboidal domain, and reformulated the Navier–Stokes equations in a non-orthogonal coordinate system using a metric tensor.

The compliant surface was modelled as a spring-backed plate, using an equation similar to (3.1) to relate the normal displacement to the pressure at the surface. However, tangential displacement was not considered in the model for the compliant wall. The authors found little modification of the near-wall turbulent structures and the drag on the surface due to wall compliance, and reported in their abstract that ‘$\ldots$ the statistical effect of the wall compliance on the turbulent channel flow is small’. A subsequent study by Kim & Choi (2014) also considered the same model, but lower spring stiffnesses and damping coefficients were considered. This study confirmed that there is very little effect on the turbulent flow for a stiff wall. For a relatively soft wall, two-dimensional downstream travelling waves are generated in the wall, which increase fluid velocity fluctuations close to the wall.

When the spring-backed wall model is used for a channel, one aspect not considered is the wall deformation due to the mean pressure gradient required to generate the flow. One implicit assumption, justified for stiff compliant walls, is that the slope of the wall is small, and there is very little change in the channel dimension over length scales comparable to the largest flow structures. A second assumption is that the wall displacement due to pressure fluctuations is much larger than the displacement due to the mean pressure difference across the length of the channel. The latter assumption is more restrictive and requires stronger justification. The wall deformation due to the shear stress tangential to the surface is also neglected, and the tangential velocity is set equal to zero at the compliant wall. The assumption here is that the wall is infinitely stiff in the tangential direction. Due to this simplification, these studies do not consider the destabilisation of the surface fluctuations due to the shear work done at the surface. Additional approximations have been made in other studies, such as the linearisation of the velocity boundary condition in Luhar, Sharma & McKeon (2015), assumption of in-plane fluctuations of the surface in Jozsa et al. (2019) or one-way coupling where the wall does not affect the fluid turbulence in Benschop et al. (2019).

A very different approach was used by Rosti & Brandt (2017) to study the flow in a channel in which one wall was made of a viscoelastic continuum. The mass and momentum equations were solved throughout the domain, and an additional order parameter which is the solid volume fraction $\phi _s$ was defined. This order parameter is $1$ in the solid and $0$ in the fluid, and the interface is identified as the location where the volume fraction passes through $\frac {1}{2}$. The volume fraction is evolved using using an advection equation. An Eulerian description is used for the solid displacement, and the evolution equation for the left Cauchy–Green tensor is obtained from the condition that its upper-convected derivative is zero. The stress is the volume-weighted sum of the solid and fluid stresses. The stress in the fluid is modelled using Newton's law of viscosity, while the neo-Hookean model is used for the solid stress. The turbulent flow simulations were conducted at a fixed Reynolds number of 2800, and the parameter $(G / \rho _f V_f^{2})$ was fixed in the range $0.5 - 2$. Rosti & Brandt (2017) found that there is substantial modification of the fluid turbulence due to the soft surface. The streamwise vorticity streaks that are observed in a turbulent flow past a rigid surface are strongly attenuated when the surface is compliant, and there are strong correlations in the spanwise direction. The maximum of the mean velocity profile is shifted away from the compliant wall as the elasticity of the wall was decreased. The logarithmic layer was modified by a shift in the origin, used for rough and permeable walls, and a change in the von Kármán constant in the logarithmic law. The shift in the velocity profile was correlated to the root mean square of the wall-normal fluctuating velocity, in a manner similar to that over rough or porous surfaces.

The simulations of Rosti & Brandt (2017) were the first to indicate that the root mean square of the velocity fluctuations and the Reynolds stress are non-zero at the fluid–solid interface. These are consistent with the experimental observations of Srinivas & Kumaran (2015, 2017b). Due to this, the wall shear stress at the compliant wall is higher than that at the rigid wall. The authors also found that the near-wall maxima in the streamwise and cross-stream fluctuating velocity and the Reynolds stress at the soft wall are higher than those near the rigid surface, indicating that the soft wall plays a role in generating turbulent fluctuations. There are large fluctuations in the wall-normal velocity at the interface between the fluid and the solid; these are comparable with the streamwise velocity fluctuations. This correlated best with the ‘wall flutter’ reported by Srinivas & Kumaran (2017b), or the ‘inviscid mode’ instability of Kumaran (1996) and Shankar & Kumaran (1999, 2000). However, in the case of wall flutter, there are visible waves on the surface of the soft solid, and the maximum of the mean velocity is closer to the compliant surface.

It should be noted that the Reynolds number used by Rosti & Brandt (2017) is much higher than that used to measure turbulent fluctuations in Srinivas & Kumaran (2015, 2017b), so a direct comparison is difficult. However, there are qualitative similarities as well as differences which require resolution. There is certainly a need for further studies over a wider parameter range in order to make a closer contact between experiments and simulations of turbulent flows.

10. Conclusions and outlook

The stability of flows in conduits with compliant walls is a relatively new field, and most of the theoretical and experimental work has been carried out in the last 25 years. The theoretical prediction of a subcritical instability in the flow past a compliant wall even in the limit of zero Reynolds number due to coupling between the wall displacement and fluid velocity (Kumaran et al. 1994; Shankar & Kumaran 2001b; Chokshi & Kumaran 2008a), discussed in § 5, and the experimental verification of this instability (Kumaran & Muralikrishnan 2000; Muralikrishnan & Kumaran 2002; Eggert & Kumar 2004), summarised in § 8.1, is a novel phenomenon uncovered in this field. At high Reynolds number, at least two different modes of instability have been identified using asymptotic analysis which are qualitatively different from the Tollmien–Schlichting modes in a rigid conduit. These are the inviscid modes with or without an internal viscous layer discussed at the end of § 6.2 and the wall mode instability (§ 6.3) which is not present in a rigid conduit (Kumaran 1998; Shankar & Kumaran 2001a, 2002; Chokshi & Kumaran 2009). As discussed in § 6.1, theorems have been derived for the possibility of an instability in the flow past compliant surfaces, similar to the Rayleigh and Fjørtoft theorems for the flow in rigid conduits for a spring-backed wall model (Yeo & Dowling 1987; Shankar & Kumaran 2000). Though the theorems are derived for the spring-backed wall model, the results appear to apply even for flows employing the viscoelastic wall model. In high-Reynolds-number experiments (§ 8.2), transition has been observed at a Reynolds number as low as 500 in a soft-walled tube (Verma & Kumaran 2012; Neelamegam & Shankar 2015) and as low as 200 in a microchannel with a soft wall (Verma & Kumaran 2013). Thus, a multi-fold reduction in the transition Reynolds number has been demonstrated for the flow conduits with compliant walls. The transition at high Reynolds number seems to be captured by an approximate local stability analysis of Verma & Kumaran (2013, 2015) (§ 8.2.2), in contrast to the flow in rigid tubes/channels where the linear stability analysis does not accurately predict transition.

The initial studies of flow past a viscoelastic wall were carried out using the linear model for solid elasticity, which is not invariant under rotation. The frame-invariant neo-Hookean model with the nonlinear strain measure (equation (3.11)) was first used by Gkanis & Kumar (2003), and this revealed an additional mode of instability due to the normal stress difference in the wall material. Consistency between the Lagrangian formulation (§ 3.3.1) and Eulerian formulation (§ 3.3.2) for the viscoelastic solid wall was shown by Chokshi & Kumaran (2008a), but the study appears to have used different constitutive relations for both formulations which coincidentally provide the same result. Subsequently, the consistency between Eulerian and Lagrangian formulations has been demonstrated by Patne et al. (2017) and Patne & Shankar (2019a) in both Cartesian and cylindrical coordinates. Though the relation between the stress and strain is still linear in the neo-Hookean and the Mooney–Rivlin models, the strain is a nonlinear function of the gradient in the displacement field. There are several intricacies in the implementation of the boundary conditions at the interface between the solid and liquid, because they have to be implemented at the displaced interface. This requires an expansion of the quantities about their values in the base state, and there are several terms in the expansion which have been missed out in previous studies. The consistent framework for the solid dynamics is non-trivial, and is one achievement of advances in the field.

One area requiring more work is the solid elasticity model that is appropriate for gels and soft tissues of relevance to experiments or practical applications. Theoretical studies so far use the neo-Hookean solid model or the more general Mooney–Rivlin solid model. In order to make quantitative comparisons for low-Reynolds-number flows where the solid strain is larger than $1$, it is necessary to have independent confirmation that the experimental stress–strain relationship is in agreement with the theoretical models used. The qualitative nature of the theoretical predictions does not depend strongly on the solid model for a Couette flow, especially when the thickness of the solid is larger than that of the fluid layer. In this case, the theoretical results are in quantitative agreement with the experimental results for the rheometer geometry shown in figure 9.

For pressure-driven flows in tubes and channels, it was shown by Gaurav (2009), Gaurav & Shankar (2010), Patne et al. (2017) and Patne & Shankar (2019a) that the results of the stability analysis strongly depend on the solid model, and even on the terms included in the interface condition between the fluid and the solid. Therefore, an accurate solid model is indispensable for an understanding of the stability in experimental systems, especially in the low-Reynolds-number limit where the solid strain in the base state is not small. Theoretical models typically assume that the solid is incompressible and the stress–strain relationship is linear in the solid, and the result for the solid deformation in the base state is unidirectional and fully developed. It is not clear if such a base state can be obtained if the stress–strain relationship is nonlinear. For quantitative comparison with linear theoretical models, it is necessary to use hyperelastic materials, which can undergo large strain, in experiments. It will also be necessary to identify materials which can undergo large strain without failing.

An important assumption in the theoretical studies conducted so far is that the solid is incompressible. In experimental systems, even if the bulk elasticity modulus is much larger than the shear modulus, it may not be accurate to assume that the solid is incompressible for pressure-driven flows. As discussed in § 3.3.6, there is a decrease in the pressure with downstream distance, and the pressure drop across the conduit is $(L/h_f)$ larger than the shear stress, where $L$ is the length and $h_f$ is the channel width/tube diameter. Thus, the pressure could be significantly larger than the shear stress for long channels/tubes, and the compressional strain could be comparable to or larger than the shear strain.

At high Reynolds number, it is shown at the end of § 3.3 that the strain in the solid in the base state is small, and the linear model for the solid elasticity is more accurate. As shown in § 8.1, the incompressibility assumption may still not be accurate in conduits where the length is much larger than the height/diameter. Therefore, a stability analysis which considers a steady and fully developed base state may not accurately predict the experimental results. Stability analyses which incorporate the deformation due to the pressure gradient have been carried out (Verma & Kumaran 2012, 2015), as indicated in § 8.2.2. Here, the deformed shape is determined from experimental measurements, and the stability analysis is carried out at different downstream locations assuming the flow is unidirectional. These assumptions are valid only if the slope of the wall is small, and the wavelength of unstable modes is smaller than the length scale for the flow development. Though results of this type of analysis are fairly accurate when compared with experiments, further progress is required. This could include more detailed models for compressible solids that can predict the deformed shape of the conduit, and global stability studies which do not employ the parallel flow approximation.

Further progress is required in identifying the mode of instability and the wavelength and frequency of the most unstable mode in experiments. There are restrictions on the conduit dimensions and the flow rate for the transition to be observed at a Reynolds number lower than the rigid-wall laminar–turbulent transition. As shown at the beginning of § 8.2, the smallest dimension of the conduit has to be about 1 mm or less, the elasticity modulus for the wall has to be 1–10 kPa and the fluid velocity has to be 1 m s$^{-1}$ or more. Due to these, the experimental geometries are small in size. The estimate for the wall displacement is of the order of tens of micrometres in typical experimental geometries, and the frequency from the theoretical predictions is 1 kHz or more. Though an increase in the root mean square displacement at transition has been detected in experiments, mapping out the displacement field of a few micrometres at high frequency in experiments is exceedingly difficult. This is an area where advances in experimental techniques are required.

Despite the unresolved issues, there is a consensus that there is one class of instability at low Reynolds number in the absence of inertia for which the transition Reynolds number scales proportional to $\varSigma$. At high Reynolds number, there are two instabilities, the inviscid modes with an internal viscous layer where the transition Reynolds number is proportional to $\varSigma ^{1/2}$ and the wall-mode instability where the transition Reynolds number is proportional to $\varSigma ^{3/4}$ (Kumaran 2003). The experimental results for the transition Reynolds number fall approximately in the range between $\varSigma ^{1/2}$ and $\varSigma ^{3/4}$, though clear scaling laws are not expected due to significant wall deformation in the experiments.

The study of turbulence in conduits with compliant walls at relatively low Reynolds number is even more recent, and the first experiments were reported within the past five years. Though the field is relatively less developed, and more verification is required, the results are intriguing and suggest that soft-wall turbulence can be considered as a new class of turbulence distinct from turbulence in the flow past rigid surfaces. Both direct numerical simulations (Rosti & Brandt 2017), discussed in § 9, and experimental studies (Srinivas & Kumaran 2015, 2017b), summarised in § 8.3, suggest that the fluid velocity fluctuations do not decrease to zero at the soft wall and that the Reynolds stress is non-zero at the wall. This indicates that the wall motion plays a role in generating turbulent velocity fluctuations.

Experimental results summarised in § 8.3 show that the turbulence energy production is a maximum at the wall. This is in contrast to turbulence in the flow past a rigid wall, where the energy production is a maximum in the near-wall region due to the lift-off and bursting of vortices. This suggests that the mechanism of turbulent energy production for soft-wall turbulence is very different from that for the flow past a rigid surface. The turbulence intensity is also very high at Reynolds numbers as small as 250–400. When scaled by the square of the mean velocity, the root mean square velocity fluctuations for the flow in a microchannel at Reynolds number 250–400 are higher than those in a rigid channel at Reynolds numbers 12 000 to 20 000. Thus, soft-wall turbulence appears to be a distinct class of turbulence where the wall dynamics appears to play a role in generating turbulence.

A distinct transition between hard-wall turbulence and soft-wall turbulence has also been observed in a channel at a Reynolds number higher than the hard-wall transition Reynolds number of about 1000 (Srinivas & Kumaran 2017b). This suggests that transition due to wall mobilisation and fluid–wall coupling occurs in a flow that is already turbulent. This transition between two distinct types of turbulence is unusual and significant. The transition Reynolds number follows the same scaling law as that for the the laminar–turbulent transition in the flow past a flexible surface. It is necessary to carry out a stability analysis of the turbulent flow past a stationary surface to determine the mechanism of transition.

There are still unresolved issues. Due to the small dimension of the channel, the spatial resolution of the velocity fields is limited. Due to this, it is has not been possible, so far, to identify coherent structures, though the mean and root mean square velocities have been measured. One important unresolved issue is the nature of the velocity field close to the wall. From the experiments, it is not clear whether a viscous sublayer exists and is not resolved, or whether a viscous sublayer is absent in soft-wall turbulence. The evolution of the displacement field in the wall has also not been resolved, though the mean square of the displacement fluctuations has been measured. A lot more work is required for a good understanding of this phenomenon.

Microfluidic mixing is an important application for the transition in compliant conduits. Microfluidic devices are currently operated in the laminar regime due to the small size of the microfluidic channels, and cross-stream mixing is due to diffusion. As estimated in § 8.5, the diffusion time across a distance of $1$ mm in a microchannel is $10^{3}$ s for small molecules in water and it could be as high as $10^{6}$ s for large and complex biomolecules. Consequently, the length of the conduits in microfluidic channels has to be of the order of tens of centimetres for complete mixing. A very large pressure drop is required for flow in lengthy channels with small cross-sectional area; this results in higher energy cost and the requirement of higher material strength of components to avoid failure. Several strategies have been proposed for enhancing mixing, but all of these have their disadvantages. Examples of passive strategies are tortuous channels (Jiang et al. 2004), sequential split-and-recombine procedures (Bessoth, de Mello & Manz 1999; Lee et al. 2006) and creation of secondary flows by wall roughness (Stroock et al. 2002). The disadvantage is that the pressure drop required is further enhanced in curved or rough channels for equal length and cross-sectional areas. In active strategies, moving parts such as pumps and stirrers are used to enhance mixing (Mensing et al. 2004). The cost involved in fabricating small intricate moving parts and the ancillary equipment required to drive these is the disadvantage of active strategies. The soft-wall instability is an optimal combination of the two, because mixing is enhanced by wall mobilisation which is due to a flow instability, and not due to external actuation. It has been shown (Kumaran & Bandaru 2016) that the mixing time is reduced by up to five orders of magnitude due to this instability, and this overcomes a significant technological barrier for on-chip sample preparation and mixing in microfluidic devices (Abbasi et al. 2019).

The transition in compliant tubes is also relevant for cardiovascular flows, where the walls of the arteries/veins are made of soft tissue. The flow is periodic, and the conduits are curved and branched. However, transition could still be of relevance, because the time scale for the growth or the instability is much smaller than the time period of the imposed flow. It has been known for some time (Aars & Solberg 1971; Stein & Sabbah 1976; Ha et al. 2018) that the extent of ‘turbulence’ in aortic flow is correlated with several pathological conditions. However, the term turbulence here often refers to intermittency and irregular fluid streamlines rather than fluid turbulence, and this is often correlated to blockage of arteries. The recent study of Saqr et al. (2020) suggests that physiological flows are turbulent even in normal large arteries at Reynolds number as low as 300, but the energy spectrum is different from the Kolmogorov spectrum. This is clearly reminiscent of the results discussed in Srinivas & Kumaran (2017b), which show that the transition to turbulence in a channel with soft walls of dimension ${\sim }1$ mm can occur at Reynolds number as low as 300. If ‘soft-wall turbulence’ due to a dynamical interaction between the fluid and the blood vessel is of relevance, it has major implications for the understanding and modelling of physiological blood flow.

For practical applications, another important finding is the transition from hard-wall to soft-wall turbulence discussed in § 8.3. This implies that fluid–wall coupling in compliant conduits could cause a transition even for flows that are already turbulent. The studies of this subject are as yet preliminary, but there is convincing evidence that there is a transition from hard-wall turbulence to a flow where the wall is mobilised and coupled to the fluid turbulence. The characteristics of soft-wall turbulence appear to be very different from those of the flow past rigid surfaces. This intriguing facet of wall-bounded turbulent flows could have far-reaching consequences in applications where the flow boundaries are sufficiently soft that there is a dynamical coupling with the fluid turbulence.