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 Reference Benjamin1960, Reference Benjamin1963; Landahl Reference Landahl1962; Carpenter & Garrad Reference Carpenter and Garrad1985, Reference Carpenter and Garrad1986; Yeo & Dowling Reference Yeo and Dowling1987), and less often for internal flows (Gajjar & Sibanda Reference Gajjar and Sibanda1996; Kumaran Reference Kumaran1996; Davies & Carpenter Reference Davies and Carpenter1997a; Larose & Grotberg Reference Larose and Grotberg1997; Shankar & Kumaran Reference Shankar and Kumaran2000; Thaokar, Shankar & Kumaran Reference Thaokar, Shankar and Kumaran2001). 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. (Reference Thaokar, Shankar and Kumaran2001) for a viscous flow and in Larose & Grotberg (Reference Larose and Grotberg1997) for the high-Reynolds-number flow in a compliant channel. Thaokar et al. (Reference Thaokar, Shankar and Kumaran2001) 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 Reference Drazin and Reid1981). An example of the derivation for a plane channel, along the lines of Yeo & Dowling (Reference Yeo and Dowling1987), is given in § 6.1. The equivalents for a pipe flow have been derived in Kumaran (Reference Kumaran1996) and Shankar & Kumaran (Reference Shankar and Kumaran2000).

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 (Reference Kumaran and Srivatsan1998) and Thaokar & Kumaran (Reference Thaokar and Kumaran2002) 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 (Reference Stewart, Waters and Jensen2009, Reference Stewart, Waters and Jensen2010), 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 (Reference Kumaran and Srivatsan1998) and Stewart et al. (Reference Stewart, Waters and Jensen2009, Reference Stewart, Waters and Jensen2010) 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 (Reference Thaokar and Kumaran2002), 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 (Reference Harden and Pleiner1994), 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 (Reference Thaokar and Kumaran2002) 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 Reference Kumaran, Fredrickson and Pincus1994; Shankar & Kumaran Reference Shankar and Kumaran2001b; Chokshi & Kumaran Reference Chokshi and Kumaran2007, Reference Chokshi and Kumaran2008a, Reference Chokshi and Kumaran2009) 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 (Reference Gaurav and Shankar2010) and Patne & Shankar (Reference Patne and Shankar2019a,Reference Patne and Shankarb) 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 Reference Kumaran1995, Reference Kumaran1996, Reference Kumaran1998; Shankar & Kumar Reference Shankar and Kumar2004; Shankar & Kumaran Reference Shankar and Kumaran1999, Reference Shankar and Kumaran2000, Reference Shankar and Kumaran2001a, Reference Shankar and Kumaran2002).

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 (Reference Malvern1969) and Holzapfel (Reference Holzapfel2000). 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 (Reference Malvern1969) and Holzapfel (Reference Holzapfel2000). 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 (Reference Chokshi and Kumaran2007, Reference Chokshi and Kumaran2008a,Reference Chokshi and Kumaranb), and is referred to as the L3 formulation in Patne, Giribabu & Shankar (Reference Patne, Giribabu and Shankar2017) and Patne & Shankar (Reference Patne and Shankar2019a). 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 (Reference Gkanis and Kumar2005), Gaurav (Reference Gaurav and Shankar2009) and Gaurav & Shankar (Reference Gaurav and Shankar2010) 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. (Reference Patne, Giribabu and Shankar2017) and Patne & Shankar (Reference Patne and Shankar2019a), 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. (Reference Patne, Giribabu and Shankar2017) and Patne & Shankar (Reference Patne and Shankar2019a).

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 Reference Gaurav and Shankar2009) 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).

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. (Reference Kumaran, Fredrickson and Pincus1994) 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 (Reference Chokshi and Kumaran2008a).

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. (Reference Kumaran, Fredrickson and Pincus1994) 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 (Reference Gkanis and Kumar2003) 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 (Reference Chokshi and Kumaran2008a) 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. Reference Kumaran, Fredrickson and Pincus1994) 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 (Reference Chokshi and Kumaran2008a) 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 (Reference Chokshi and Kumaran2008a) 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 (Reference Chokshi and Kumaran2008a) 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 (Reference Shankar and Kumaran2001b), 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 (Reference Chokshi and Kumaran2008a), 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. (Reference Thaokar, Shankar and Kumaran2001) 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 (Reference Thaokar and Kumaran2002), 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 Reference Kumaran1995) 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 (Reference Gkanis and Kumar2005), 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 (Reference Gaurav and Shankar2010) 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 (Reference Gaurav and Shankar2009) 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 (Reference Gkanis and Kumar2005) and Gaurav & Shankar (Reference Gaurav and Shankar2010) 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 (Reference Gkanis and Kumar2005), but was included in Gaurav & Shankar (Reference Gaurav and Shankar2010). On the solid side, Gkanis & Kumar (Reference Gkanis and Kumar2005) 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 (Reference Gaurav and Shankar2010) because the boundary conditions were applied at the undeformed interface in the Lagrangian formulation. Neither Gkanis & Kumar (Reference Gkanis and Kumar2005) nor Gaurav & Shankar (Reference Gaurav and Shankar2010) 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 (Reference Gaurav and Shankar2010) and Gaurav (Reference Gaurav and Shankar2009) 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. (Reference Patne, Giribabu and Shankar2017) and Patne & Shankar (Reference Patne and Shankar2019a), 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.1. Theorems for inviscid flow

Equivalents of the classical theorems for inviscid flows (Drazin & Reid Reference Drazin and Reid1981) 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 (Reference Drazin and Reid1981) and the adaptation for a channel with a compliant wall is given in Yeo & Dowling (Reference Yeo and Dowling1987). The derivation of the classical theorems for a rigid tube was first carried out by Howard & Gupta (Reference Howard and Gupta1962) and Maslowe (Reference Maslowe1974), and the calculation for a tube with a compliant wall is given in Kumaran (Reference Kumaran1996) and Shankar & Kumaran (Reference Shankar and Kumaran2000).

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 Reference Yeo and Dowling1987) 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 (Reference Shankar and Kumaran2000) 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 (Reference Reynolds1883).

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 (Reference Drazin and Reid1981). The matched asymptotic analysis was first carried out in the classic works of Tollmien (Reference Tollmien1929) and Lin (Reference Lin1945), 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 Reference Corcos and Sellars1959; Gill Reference Gill1965) 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 (Reference Kumaran1998) 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 Reference Shankar and Kumaran2001a, Reference Shankar and Kumaran2002). 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 Reference Kumaran1998; Shankar & Kumaran Reference Shankar and Kumaran2001a, Reference Shankar and Kumaran2002). 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 (Reference Rotenberry1992), Rotenberry & Saffman (Reference Rotenberry and Saffman1990), Davies & Carpenter (Reference Davies and Carpenter1997a,Reference Davies and Carpenterb) and Larose & Grotberg (Reference Larose and Grotberg1997). In these studies, with the exception of Larose & Grotberg (Reference Larose and Grotberg1997), 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 (Reference Rotenberry1992), Rotenberry & Saffman (Reference Rotenberry and Saffman1990) and Davies & Carpenter (Reference Davies and Carpenter1997b), 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 (Reference Rotenberry and Saffman1990) 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 (Reference Davies and Carpenter1997a). 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 Reference Carpenter and Garrad1985, Reference Carpenter and Garrad1986). 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 (Reference Huang1998) 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 (Reference Gajjar and Sibanda1996), 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 (Reference Larose and Grotberg1997), 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 (Reference Larose and Grotberg1997) 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. (Reference Stewart, Waters and Jensen2010). 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. (Reference Stewart, Waters and Jensen2009) 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 (Reference Kumaran and Srivatsan1998) 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 (Reference Chokshi and Kumaran2008b), 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 (Reference Pierce1992) 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 (Reference Kumaran1996) 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 (Reference Shankar and Kumaran1999) 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 (Reference Shankar and Kumaran2000), 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 (Reference Kumaran1996) and Shankar & Kumaran (Reference Shankar and Kumaran1999, Reference Shankar and Kumaran2000) 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 (Reference Gaurav and Shankar2009) and Gaurav & Shankar (Reference Gaurav and Shankar2010). 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. (Reference Patne, Giribabu and Shankar2017) and Patne & Shankar (Reference Patne and Shankar2019a), 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.

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 (Reference Srinivas and Kumaran2017b) 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 (Reference Kumaran and Muralikrishnan2000) and Muralikrishnan & Kumaran (Reference Muralikrishnan and Kumaran2002), followed by Eggert & Kumar (Reference Eggert and Kumar2004) on the dynamics after transition and Neelamegam, Giribabu & Shankar (Reference Neelamegam, Giribabu and Shankar2014) 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 (Reference Kumaran and Muralikrishnan2000).

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 (Reference Kumaran and Muralikrishnan2000) and Neelamegam et al. (Reference Neelamegam, Giribabu and Shankar2014). A comparison of the experimental result of Kumaran & Muralikrishnan (Reference Kumaran and Muralikrishnan2000) with the theoretical prediction of Kumaran et al. (Reference Kumaran, Fredrickson and Pincus1994) 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. (Reference Neelamegam, Giribabu and Shankar2014) 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 (Reference Kumaran and Muralikrishnan2000) and Muralikrishnan & Kumaran (Reference Muralikrishnan and Kumaran2002) established that there is an instability at a critical strain rate, the nature of the flow after transition was not explored. Eggert & Kumar (Reference Eggert and Kumar2004) 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 (Reference Eggert and Kumar2004) 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 (Reference Eggert and Kumar2004) is in agreement with the predictions of the weakly nonlinear stability analysis of Shankar & Kumaran (Reference Shankar and Kumaran2001b) and Chokshi & Kumaran (Reference Chokshi and Kumaran2008a).

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 (Reference Eggert and Kumar2004).

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 (Reference Patel and Head1969), 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 (Reference Darbyshire and Mullin1995), 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 (Reference Verma and Kumaran2012, Reference Verma and Kumaran2013) and Neelamegam & Shankar (Reference Neelamegam and Shankar2015), 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 (Reference Verma and Kumaran2012), Neelamegam & Shankar (Reference Neelamegam and Shankar2015) and Chandra, Shankar & Das (Reference Chandra, Shankar and Das2019) 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 (Reference Verma and Kumaran2012), about 1.65 mm in Neelamegam & Shankar (Reference Neelamegam and Shankar2015) and about 0.4 mm in Chandra et al. (Reference Chandra, Shankar and Das2019). 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 (Reference Verma and Kumaran2012) and Neelamegam & Shankar (Reference Neelamegam and Shankar2015) 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 (Reference Verma and Kumaran2012) and Neelamegam & Shankar (Reference Neelamegam and Shankar2015), 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 (Reference Verma and Kumaran2012) 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. (Reference Chandra, Shankar and Das2019) 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. (Reference Chandra, Shankar and Das2019) 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 (Reference Verma and Kumaran2012) and Neelamegam & Shankar (Reference Neelamegam and Shankar2015) and as low as $250$ in the study of Chandra et al. (Reference Chandra, Shankar and Das2019) 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 (Reference Verma and Kumaran2012) 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 (Reference Verma and Kumaran2012).

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 (Reference Verma and Kumaran2012).

In order to rule out the effect of the change in shape of the tube on the transition, Verma & Kumaran (Reference Verma and Kumaran2012) 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 Reference Reynolds1883), 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 (Reference Verma and Kumaran2012) and Neelamegam & Shankar (Reference Neelamegam and Shankar2015). 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 (Reference Verma and Kumaran2012) 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 (Reference Verma and Kumaran2012), and this provided the scaling law ${Re}_t \propto \varSigma ^{5/8}$ for the wall mode instability. Neelamegam & Shankar (Reference Neelamegam and Shankar2015) 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 (Reference Neelamegam and Shankar2015) 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. (Reference Chandra, Shankar and Das2019) 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 (Reference Verma and Kumaran2012).

Figure 14. The transition Reynolds number as a function of the parameter $\varSigma$ from Verma & Kumaran (Reference Verma and Kumaran2012) (black), Neelamegam & Shankar (Reference Neelamegam and Shankar2015) (blue) and Chandra et al. (Reference Chandra, Shankar and Das2019) (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 (Reference Verma and Kumaran2012).

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 (Reference Verma and Kumaran2015). 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 (Reference Verma and Kumaran2015) 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 (Reference Verma and Kumaran2015) 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 (Reference Verma and Kumaran2013). 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 (Reference Verma and Kumaran2013).

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 (Reference Verma and Kumaran2013).

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 (Reference Verma and Kumaran2013).

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 (Reference Verma and Kumaran2013). Verma & Kumaran (Reference Verma and Kumaran2013) 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 (Reference Verma and Kumaran2013). (b) The transition Reynolds number from Verma & Kumaran (Reference Verma and Kumaran2013) (blue), Srinivas & Kumaran (Reference Srinivas and Kumaran2017b) (black) and Kumaran & Bandaru (Reference Kumaran and Bandaru2016) (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 (Reference Verma and Kumaran2013).

The transition Reynolds numbers from the study of Verma & Kumaran (Reference Verma and Kumaran2013) are plotted along with those of Srinivas & Kumaran (Reference Srinivas and Kumaran2017b) 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 (Reference Kumaran and Bandaru2016), 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 (Reference Srinivas and Kumaran2017b) 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 (Reference Srinivas and Kumaran2015) and Srinivas & Kumaran (Reference Srinivas and Kumaran2017b), 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 (Reference Srinivas and Kumaran2015) was carried out for the same set-up as Verma & Kumaran (Reference Verma and Kumaran2013). 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 (Reference Srinivas and Kumaran2015).

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 (Reference Srinivas and Kumaran2015), 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 (Reference Srinivas and Kumaran2015).

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 (Reference Srinivas and Kumaran2015) 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 (Reference Srinivas and Kumaran2017b) 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 (Reference Srinivas and Kumaran2017b). 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 (Reference Srinivas and Kumaran2017b) 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 (Reference Srinivas and Kumaran2017b).

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 (Reference Srinivas and Kumaran2017b).

The study of Srinivas & Kumaran (Reference Srinivas and Kumaran2017b) 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 (Reference Srinivas and Kumaran2017b) 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 (Reference Neelamegam, Shankar and Das2013) in the context of the ‘purely elastic’ instability due to the hoop stress in a rotating polymeric fluid (Shaqfeh Reference Shaqfeh1996; Groisman & Steinberg Reference Groisman and Steinberg2004). 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. (Reference Neelamegam, Shankar and Das2013) 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 (Reference Srinivas and Kumaran2017a) for the flow through a PDMS microchannel shown in figure 15 and by Chandra et al. (Reference Chandra, Shankar and Das2019) 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 (Reference Srinivas and Kumaran2017a), and in the range 30–50 kPa in Chandra et al. (Reference Chandra, Shankar and Das2019). 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 (Reference Srinivas and Kumaran2017a), 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. (Reference Chandra, Shankar and Das2019) 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 Reference Doi and Edwards1988), first normal stress difference (Larson, Shaqfeh & Muller Reference Larson, Shaqfeh and Muller1990) and capillary break-up (CaBeR) technique for the extensional relaxation time (Dinic et al. Reference Dinic, Zhang, Jimenez and Sharma2015). Zell et al. (Reference Zell, Gier, Rafai and Wagner2010) found that there are significant differences in the results of different measurements. Srinivas & Kumaran (Reference Srinivas and Kumaran2017a) 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. (Reference Chandra, Shankar and Das2019) 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 (Reference Srinivas and Kumaran2017a) 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. (Reference Chandra, Shankar and Das2019), 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. (Reference Chandra, Shankar and Das2019) 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 (Reference Srinivas and Kumaran2017a) 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 Reference de Gennes1976, Reference de Gennes1979). 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 (Reference Srinivas and Kumaran2017a) were entirely in the dilute regime, whereas the experiments of Chandra et al. (Reference Chandra, Shankar and Das2019) extended over the dilute and semi-dilute regimes.

The experimental set-up of Srinivas & Kumaran (Reference Srinivas and Kumaran2017a) was the same as that of Verma & Kumaran (Reference Verma and Kumaran2013) 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 (Reference Srinivas and Kumaran2017a) 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. (Reference Chokshi, Bhade and Kumaran2015) for the effect of viscoelasticity on the wall mode instability in a channel.

The second important observation in Srinivas & Kumaran (Reference Srinivas and Kumaran2017a) 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. (Reference Chandra, Shankar and Das2019) 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 (Reference Verma and Kumaran2012). Chandra et al. (Reference Chandra, Shankar and Das2019) concluded that the destabilisation is due the wall mode instability modified by fluid viscoelasticity, as reported by Srinivas & Kumaran (Reference Srinivas and Kumaran2017a). For higher polymer concentration, Chandra et al. (Reference Chandra, Shankar and Das2019) found a very different cause of destabilisation, which is the elasto-inertial instability (Zakin et al. Reference Zakin, Ni, Hansen and Reischman1977; Samanta et al. Reference Samanta, Dubief, Holzner, Schafer, Morozov, Wagner and Hof2013) 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 Reference Chandra, Shankar and Das2018), 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. (Reference Chandra, Shankar and Das2019) 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 (Reference Shrivastava, Cussler and Kumar2008) 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 (Reference Verma and Kumaran2013) 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 (Reference Verma and Kumaran2012) 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 (Reference Kumaran and Bandaru2016) 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 (Reference Kumaran and Bandaru2016) 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. (Reference Abbasi, Chowdhury, Subramaniam, Jain, Muthe, Sheikh, Banerjee and Kumaran2019), the utility of this technology has recently been demonstrated for a point-of-care device.