Stability of fluid flows coupled by a deformable solid layer

Abstract Fluid flows coupled by a deformable solid layer (DS) may be found in both natural and industrial settings. To analyse the stability of such systems, we consider plane Couette flows coupled by a DS. The shear stress exerted by the flowing fluids on the DS leads to its deformation, coupling the motion in each fluid. We demonstrate that, as a consequence of the coupling and irrespective of the thickness of the DS, a viscous instability exists, which can lead to absolute instability at sufficiently high dimensionless speed of the lower plate. The predicted instability apparently arises due to the energy exchange between the fluids and DS via the tangential velocities at the interface. A linear elasticity model for the DS is employed, removing the non-physical growth rate predicted by previous models. A spatio-temporal analysis further reveals the existence of an absolute instability. Insight on this instability could be potentially utilised in manipulating mixing in microfluidic systems and membrane separation as well as for understanding morphologies evolving during membrane formation via interfacial polymerisation.


Introduction
Viscous flows coupled by a deformable solid layer (DS) may be encountered in both natural and engineered processes. In biological systems, such flows occur at the micro-scale, past deformable lipid membranes forming the cell wall and organelles (Fricke 1925;Hochmuth, Mohandas & Blackshear 1973;Thaokar & Kumaran 2002), as well as at a macroscopic level, such as the heart, where the interventricular and interatrial septums act as a DS (Rosenquist et al. 1979;Marcomichelakis et al. 1983;Pislaru et al. 2014). Engineered systems where such flows are present include membrane separation (Strathmann 2011) and interfacial polymerisation (Ukrainsky & Ramon 2018) processes.
The experiments of Kumaran & Muralikrishnan (2000), Srinivas & Kumaran (2017), Neelamegam & Shankar (2015) and Verma & Kumaran (2012 showed that the presence of a deformable wall in a micro-channel can lead to transitional flows at a much lower Reynolds number compared to rigid-walled channels and, presumably, lead to improved mixing. However, this requires a wall much thicker than the channel dimensions. Furthermore, manipulation of the hydrodynamic instabilities in such a system requires alteration of the fluid and/or deformable wall properties (Verma & Kumaran 2012. Thus, during the process, mixing in such a system cannot be actively controlled. Here, we propose that such limitations can be overcome by using a thin, elastic membrane as a DS along with an adjacent, second fluid, as illustrated in figure 1. At any instant during a process, the speed of fluid 2 can be used to actively control and trigger hydrodynamic instabilities which in turn can be used to control the mixing in the coupled system. The present theoretical analysis can help in understanding the manipulation of the hydrodynamic instabilities in such systems. Such a system is also representative of some membrane-based separation processes using hollow fibres (e.g. dialysis), essentially a tube whose wall is a permeable, elastic polymer (Strathmann 2011). This membrane separates two liquid streams flowing at either side. As with the microfluidic system, flow-induced wall vibrations can enhance mixing and intensify the process. Finally, another system of interest is the stability of thin films generated during interfacial polymerisation, commonly used to fabricate desalination membranes (Raaijmakers & Benes 2016). Here, the ultra-thin polymeric film is formed at a liquid-liquid interface, and its ultimate morphology has been linked to several possible mechanisms, one of which may be of a hydrodynamic origin (Ukrainsky & Ramon 2018); the existence of such a mechanism is also part of the motivation to study the stability of this system.

Flow past a deformable solid
Before proceeding with the analysis, we begin with an overview of previous studies on related systems. The first study concerning the stability of fluids coupled by a DS was done by Harden & Pleiner (1994), where the DS was treated as an infinitesimally thin elastic sheet. Their study was aimed at explaining experimental studies on dynamic light scattering from insoluble polymeric monolayers at a liquid-liquid interface. Accordingly, the study was focused on hydrodynamic modes induced by thermal fluctuations. As such, their model neglected the governing equations for the solid film. The elastic character of the film affected the liquid-liquid interface only through the tangential stress continuity condition. Other interfacial conditions were those normally applied for a liquid-liquid interface. Additionally, both fluids were static, hence the effect of shear on the membrane and hydrodynamic modes was not studied, in contrast to our present motivation. Kumaran & Srivatsan (1998) added fluid shear and considered the stability of plane Couette flows coupled by a (still) infinitesimally thin DS. Similar to Harden & Pleiner (1994), in their study, the DS was sandwiched between two fluids as shown in figure 1. However, Kumaran & Srivatsan (1998) assumed only normal displacement of the membrane and neglected possible tangential displacement. As a result, a constitutive equation for the DS became unnecessary. Their analysis predicted that in the absence of inertia, the system is stable and the transition Reynolds number is proportional to the dimensionless tension in the thin elastic sheet.
The lack of a constitutive equation in the analysis of Kumaran & Srivatsan (1998) was then relaxed by Thaokar & Kumaran (2002), who used a constitutive relation, originally proposed by Harden & Pleiner (1994), for the DS. The flow geometry in Thaokar & Kumaran (2002) was similar to the one considered in figure 1. With the added constitutive relation, Thaokar & Kumaran (2002) predicted the existence of an instability even in the absence of inertia, in contradistinction with Kumaran & Srivatsan (1998). Thaokar & Kumaran (2002) also performed a weakly nonlinear stability analysis and predicted that perturbations in the limit of vanishing wavenumber are supercritically stable. However, similar to Harden & Pleiner (1994) and Kumaran & Srivatsan (1998), Thaokar & Kumaran (2002) assumed an infinitesimally thin elastic sheet. As shown later in § 3, their constitutive equation and/or the assumption of infinitesimally thin membrane leads to a seemingly non-physical growth rate in some parametric regimes. This prompted the use of a linear elasticity model and a finite-thickness DS in the present work, consequences of which are discussed in § 3.

Spatio-temporal stability analysis
The previous studies by Harden & Pleiner (1994), Kumaran & Srivatsan (1998) and Thaokar & Kumaran (2002) only considered the temporal evolution of disturbances in the flows coupled by a DS. However, a temporal stability analysis does not provide information about the growth of the disturbances in space or in both space and time simultaneously. This requires a spatio-temporal stability analysis, which may be further classified into either absolute or convective instabilities. A convective instability implies that, given sufficient time, the disturbances will decay at any point in space, while an absolute instability leads to growth of disturbances at any point in space. Further description of both absolute and convective instabilities, as well as an outline of the methodology used to determine their existence, is presented in § 3.2.
An absolute instability has been experimentally observed in several fluid flows, for example, Guillot et al. (2007) and Utada et al. (2008) observed that an absolute instability leads to jet breakup. The experiments of Lingwood (1995) on the boundary layer on a rotating disk showed that the laminar to turbulent flow transition could be a consequence of an absolute instability. In a rotating Hagen-Poiseuille flow, Shrestha et al. (2013) observed the presence of a wavy pattern in the whole tube, presumably due to an absolute instability. Thus, if an absolute instability exists in a given flow geometry, it can bring about dramatic changes. This feature could be exploited in promoting mixing. The length of the system may not be sufficiently large for convective instabilities to grow within the domain of the device. However, if absolute instabilities are present, then the flow is destabilised over the entire domain of the channel, which can enhance mixing, providing further motivation and potentially important practical consequences for probing the existence of absolute instabilities in fluid flows coupled by a DS.
In the present study, linear elasticity is used to model the DS due to the following reasons. For high dimensionless thickness of the DS, a, the study of Gkanis & Kumar (2003) on the plane Couette flow past a neo-Hookean solid shows that the linear elastic and neo-Hookean models are in excellent agreement. However, at low a, for the plane Couette flow past a deformable solid, the critical dimensionless speed of the plates, Γ , which is proportional to the dimensionless speed of the moving plate, is sufficiently high such that deviations appear in results obtained by linear elasticity, compared with a neo-Hookean model. However, this is not an issue in the present analysis since the coupling instability predicted here for a < 1 exists at low values of Γ (< 0.1), for which linear elasticity and more advanced models such as the neo-Hookean model are in good agreement (Gkanis & Kumar 2003;Shankar & Kumar 2004;Gaurav Shankar 2009). Additionally, the linear elasticity model is mathematically and numerically simpler compared to the neo-Hookean model. As shown later in § 3, the linear elastic model removes the unphysical growth rate predicted by the model of Thaokar & Kumaran (2002). Therefore, linear elasticity is deemed sufficient to define the dynamics of the DS and hence to study the linear stability of the present system.
The rest of the paper is arranged as follows. The base-state variables and perturbation governing equations are derived in § 2. The temporal stability results and their discussion is given in § 3.1. The existence of the absolute instability and its dependence on the various parameters is discussed in § 3.2. The salient conclusions of the present study are summarised in § 4.

Problem formulation
The system under consideration consists of an incompressible, isotropic and homogeneous elastic solid with shear modulus G, sandwiched between two incompressible fluids, as shown in figure 1. The fluids, marked as 1 and 2, have viscosities μ 1 and μ 2 , respectively and density ρ. The densities of the two fluids and the DS are assumed to be equal, which is a reasonable approximation for the industrial processes and biological systems motivating this problem (Hochmuth et al. 1973;Strathmann 2011;Shrestha et al. 2013;Pislaru et al. 2014). The lengths in the present problem are scaled by the thickness of fluid 1, R. In the dimensionless coordinates, fluid 1, DS and fluid 2 are present in the domains [0, 1], [1, 1 + a] and [1 + a, 1 + a + b], respectively. The length and height of the system are assumed to infinitely extend in the x and z-directions. The plates at y = 0 and y = 1 + a + b are moving at dimensionless speeds 1 and V r , respectively, where velocities are scaled by the speed of plate at y = 0, V 1 . The parameter V r = V 2 /V 1 is the ratio of the plate velocities where V 2 is the dimensional velocity of the plate at y = 1 + a + b.
The DS is assumed to be pre-stretched due to the tethering of its ends, in line with practical applications. The tension induced in the DS as a consequence of the tethering is assumed to be much higher than the base state stresses exerted by the flowing fluids. If the DS is un-tethered, then the present analysis is applicable if the interfacial tension is sufficiently large. Otherwise, the fluid shear stress would induce a pressure gradient in the solid that would, in turn, be transmitted to the fluid. Under such conditions, a fully developed Couette flow would no longer be a correct description of the velocity field. A force balance at the fluid-DS interface shows that, for the development of the plane Couette flow for the un-tethered DS, the interfacial tensions must be higher than (Γ /ak 3 )((V r μ r /b) + 1), where k is the wavenumber and Γ is the dimensionless speed of the plate at y = 0. The parameter Γ = μ 1 V 1 /(GR) can be interpreted as the ratio of the viscous stresses in fluid 1 to the elastic stresses within the DS. Another interpretation of Γ is the ratio of the relaxation time scale of the DS (∼ μ 1 /G) to the convection time scale (∼ R/V 1 ). The latter interpretation indicates that Γ is a measure of the relative elasticity of the DS. Furthermore, it is assumed that the dimensional widths of each fluid domain, namely R and bR, satisfy the condition R L and bR L, where L is the length of the DS. This assumption is necessary for using the normal mode analysis employed in the present study. However, if the widths of the fluids and length of the DS are of comparable magnitude, then a fully global stability analysis is required (Theofilis 2011;Tsigklifis & Lucey 2017).
We consider two fluids with a velocity field v (i) = (v (i) x , v (i) y ), where i = 1, 2 represent fluids 1 and 2, respectively. Scaling the time, t, by R/V 1 , velocities by V 1 , the dimensionless continuity equation, for an incompressible flow, becomes where the gradient operator is ∇ = e x (∂/∂ x) + e y (∂/∂ y) with e x and e y the unit vectors in the x and y-directions, respectively. The dimensionless Navier-Stokes equation for the fluids, with μ 1 V 1 /R used for scaling the pressure, is, where Re = ρV 1 R/μ 1 is the Reynolds number. The dimensionless viscosity is μ (i) = 1 for fluid 1 and μ (i) = μ r for fluid 2. Assuming a steady-state and fully developed flow, the base-state velocities of the fluids arev (1) where an overbar indicates a base-state quantity. The subsequent linear stability analysis is performed for this base state. For the DS, a linear elastic model is used. Scaling the DS displacements by R and pressure by μ 1 V 1 /R, the dimensionless incompressibility and linear momentum equations for the DS are Here, u = (u x , u y ) is the displacement field in the DS with u x and u y the x and y displacement components, respectively, and p m is the pressure field in the DS.

Linearised perturbation equations
For the linear stability analysis, two-dimensional disturbances are a reasonable assumption since for the flows past a deformable surface, two-dimensional disturbances have been shown to be more unstable than three-dimensional disturbances (Patne & Shankar 2018). For the purpose of the linear stability analysis, dynamical quantities such as velocities, displacements and pressures are decomposed into the base state and perturbed state, as . Here, f (x, t) is any dynamic quantity and a prime signifies the small perturbation quantity. In the linearised governing equations, the normal modes of the following form are then substituted where k = k r + ik i is the complex wavenumber andf ( y) is the eigenfunction of f (x, t).
The other parameter, ω = ω r + iω i is the complex frequency, which characterises the temporal frequency and growth of the disturbances. Therefore, the flow is considered to be temporally unstable if at least one eigenvalue satisfies the condition ω i > 0. For the existence of the spatio-temporal instability, as explained in § 3.2, the cusp point coordinates in the ω plane must possess ω i > 0. After substitution of the normal modes, the linearised governing equations for the fluid become (2.13) The above equations are to be solved using the following boundary conditions. At y = 0 and y = 1 + a + b, the assumption of no slip and impermeability of the plates gives v (1) At y = 1, oscillations of the fluid-DS interface will be induced due to the perturbations. The perturbed interface is y = u y (x, t)| y=1 so that the linearised normal to the interface is n (1) = −e x (∂u y (x, t)| y=1 /∂ x) + e y , pointing in the positive y-direction i.e. from the fluid 1 into the DS, and the tangent is t (1) = e y (∂u y (x, t)| y=1 /∂ x) + e x . There will be four boundary conditions at y = 1, two of which are continuity of the tangential and normal velocities, while the other two represent continuity of the tangential and normal stresses. These boundary conditions, after substitution of the normal modes, becomẽ where all quantities are evaluated at y = 1. The parameter Σ 1 is the dimensionless interfacial tension between fluid 1 and the DS. Similarly, for the interface at y = 1 + a, the equation of the surface is y = −u y (x, t)| y=1+a . Thus, the linearised normal and tangent to the fluid-DS interface are n (2) = e x (∂u y (x, t)| y=1+a /∂ x) + e y and t (2) = −e y (∂u y (x, t)| y=1+a /∂ x) + e x , respectively. To maintain consistency with the direction of n (1) , n (2) is directed from the fluid 2 into the DS. The boundary conditions where all quantities are evaluated at y = 1 + a. Finally, Σ 2 is the dimensionless interfacial tension between fluid 2 and the DS. The term Dv (i) xũ y in the tangential velocity continuity conditions at the interface, (2.16) and (2.20), leads to the coupling between the base state and perturbations, which in turn leads to the energy exchange between the DS and the fluids. This energy exchange can eventually lead to the instabilities predicted in the present work. However, we also note that, at higher order, the interaction between the normal interfacial displacement and the pressure fluctuations represents another mode of energy transfer. To determine the stability of the system, the above equations and boundary conditions are solved for the eigenvalue ω using the numerical methodology outlined in the next subsection.

Numerical methodology
To carry out the linear stability analysis of the present problem for arbitrary Re, the pseudo-spectral method is employed, in which the eigenfunctions are expanded in terms of Chebyshev polynomials, asf (2.24) wheref ( y), m, N, T m ( y) and a m are, respectively, the eigenfunctions, the number of the Chebyshev polynomial, the highest degree of the polynomial in the series expansion or the number of the collocation points, the mth Chebyshev polynomial and the coefficient of the latter in the expansion. The series expansions are evaluated at N collocation points to determine the series' coefficients a m and/or to obtain the eigenvalues ω. The eigenvalue problem may then be written in the form where A, B and C are the discretised matrices and e is the eigenvector formed by the eigenfunctions of the dynamic quantities such as the velocity, displacements, and the fluid and DS pressures. Next, the polyeig MATLAB routine is used to solve the eigenvalue problem (2.8)-(2.23). To filter out the spurious modes from the numerically computed spectrum of the problem, the latter is obtained for N and N + 2 collocation points and the eigenvalues are compared with a specified tolerance, e.g. 10 −4 . The genuine eigenvalues are verified by increasing the number of collocation points by 20 and monitoring the variation of the obtained eigenvalues. If the eigenvalue does not change up to the sixth significant digit, then the same number of collocation points are used to predict the critical parameters. In the present work, N = 40 is found to be sufficient to determine the most unstable converged eigenvalue, within the parameter range studied.
Other parameters such as densities and viscosity of the fluids are ∼10 3 kg m −3 and ∼10 −3 − 10 −1 Pa s, respectively.
In the creeping-flow limit (Re = 0), the governing equations for the fluids and DS can be analytically solved to give eight integration constants. The obtained eigenfunctions for velocities, displacements, as well as the fluid and DS pressures are then substituted into the boundary conditions, resulting in eight equations with eight unknown variables. In matrix form, the resulting eigenvalue problem becomes where M is an 8 × 8 matrix containing coefficients of the integration constants and α is a vector containing eight integration constants. The dispersion relation, ω = ω(k, a, b, Γ, μ r , V r , Σ 1 , Σ 2 ) is then obtained by taking the determinant of the matrix M obtained in the previous step, thereby determining the system stability. In the present work, most of the results are for Re = 0 and the analytical method does not produce spurious eigenvalues, thus the analytical method is preferred for obtaining results at Re = 0.

Results and discussion
3.1. Temporal stability (k i = 0) Typical ranges of dimensional parameters involved in flows coupled by a DS, motivating this study, are shown in table 1. According to these, the Reynolds number is estimated to be Re ∼ 10 −4 − 1. Owing to such small value of the Reynolds number, Thaokar & Kumaran (2002) assumed the creeping-flow limit in their analysis. Here, we begin with the vanishing Re limit, and separately present the temporal stability results for finite Re in § 3.1.2.
Kumaran (2002) may not have been able to observe the inconsistent results predicted by their model since a counter-flow arrangement was not studied in their work. The unbounded growth rate shown in figure 2(a) implies some inconsistency in the constitutive equation or assumption of the infinitesimally thin DS. Consequently, in the present work, a linear elastic constitutive equation is employed and a finite thickness of the DS is considered. Figure 2(b) illustrates the removal of the non-physical growth rate for the most unstable (or least stable) mode by the present model, suggesting that considering a finite thickness of the DS and a linearly elastic constitutive equation, resolves the discrepancy. From figure 2(b), the growth rate peaks at k ∼ O(0.1), indicating that the dominant disturbances have a wavelength larger than the width of fluid 1. In the creeping-flow limit, the present analysis predicts four eigenvalues, shown in table 2. The instability predicted herein is due to the first eigenvalue in the table, which is a downstream travelling mode that remains downstream within the parameter range studied ('downstream mode' and 'upstream mode' here implies disturbances having positive and negative ω r , respectively). Hence, this mode is the focus of the forthcoming discussion. The remaining three modes are always stable.   x .
The form of the perturbations for the unstable eigenmode, at marginally stable parameter values, is shown in figure 3, normalised by the maximum absolute value of the corresponding eigenfunctions. In figure 3(a), at y = 0, v (1) x vanishes, as required by the boundary conditions at the moving plate near fluid 1. This is likewise observed for v (2) x , where the disturbances must vanish at y = 1 + a + b, which corresponds to y = 1 in figure 3(c). Interestingly, the absolute values achieved by the normalised perturbations v (1) x and v (2) x at the fluids-DS interface are equal. Consequently, the DS exhibits anti-symmetric oscillations (figure 3b). Thus, the unstable mode predicted in the present work is anti-symmetric in nature.
The oscillations of the DS, in the perturbed state, are shown in figure 3(b). Figure 4 shows that the second mode is symmetric with respect to the horizontal perturbation displacement of the DS, while the oscillations of the DS for the third and fourth modes are almost identical except that the fourth mode exhibits tilted oscillations. Interestingly, the oscillations for both the third and fourth modes exhibit slight variations in x and y.
In a typical experiment aimed at detecting hydrodynamic instabilities, the growth of disturbances is observed as a function of the slowly increasing velocity. In the present study, Γ is the dimensionless speed of the plate and also a measure of the shear stress applied by the flow on the DS. Therefore, we study the variation of the critical value of Γ (Γ c ), as affected by different model parameters, as an indicator of system stability. To estimate the critical parameters, neutral stability curves are obtained, as illustrated   figure 5, the neutral stability curve is flatter compared with the one for a = 0.125. Physically, this implies that at the critical speed of the plate, for a = 10, disturbances with a broader wavenumber range will become unstable. Neutral stability curves also show that with a decrease in a, Γ c decreases, demonstrating a destabilising effect when the thickness of the DS is decreased.
The critical wavenumber, k c , and Γ c , obtained from neutral stability curves calculated for various values of a, are plotted in figure 6. The figure shows a non-monotonic variation in both k c and Γ c with a decrease in a. Evidently, a decreasing thickness of the DS leads to the lowering of the instability threshold and the wavelength of the unstable disturbances increases. Figure 6 also shows the existence of an unconditionally unstable region for a < 0.12. This region is characterised by the absence of well-defined critical parameters.
For higher values of a, the curve for Γ c becomes akin to the problem of the plane Couette flow past a plate lined by a DS (Kumaran, Fredrickson & Pincus 1994), and Γ c decreases with increasing a. The plane Couette flow past a plate lined by a DS exhibits a 'viscous mode' of instability in the creeping-flow limit (Kumaran et al. 1994;Gkanis & Kumar 2003, 2005. We recall that the present system differs from the one studied by Kumaran et al. (1994) through the existence of the second fluid. Considering the variation of Γ c with a shown in figure 6, the instability is similar to the viscous mode for higher a. However, as a is decreased, the viscous mode of the instability is affected by the coupling with   xũy terms in the tangential velocity continuity equations (2.16) and (2.20). The eigenvalues are obtained at k = 0.05, a = 10 −4 , Γ = 10 −5 , b = 2, μ r = 0.5, V r = 0, Σ 1 = 1 and Σ 2 = 1. The table shows that Dv (i) xũy terms are responsible for the unconditionally unstable region, as these terms couple the base state with the perturbations. the second fluid, which leads to a rapid decrease in the critical parameters, as shown in figure 6. To conclude, the flow geometry considered here exhibits an instability similar to the viscous instability at high a; however, at low a, the coupling between the fluids via the DS modifies the viscous instability and results in an unconditionally unstable region.
Earlier studies in the creeping-flow limit, by Kumaran et al. (1994) and Gkanis & Kumar (2003), concluded that the viscous instability was due to the amplification of the disturbances at the interface between the fluid and DS. These are driven by the interaction term between the base-state velocity gradient and the perturbation in the vertical displacement of the DS, in the tangential velocity continuity equations (2.16) and (2.20). This interaction term couples the base state and perturbations, facilitating the transfer of energy from the base state to the perturbations. Further analysis reveals that the interaction terms in the tangential velocity continuity equations (2.16) and (2.20), are also the reason for the unconditionally unstable region predicted in figure 6, as shown in table 3. The existence of the unconditionally unstable region implies that as long as there is a non-zero deformation (i.e. Γ / = 0) in the DS due to the shear applied by the flow past the interface, then the coupled system is always linearly unstable provided that a and μ r are sufficiently small. Furthermore, as illustrated in table 3, this is due to the energy exchange occurring at the fluids-DS interface through Dv (i) xũ y terms in the tangential velocity continuity equations (2.16) and (2.20). We further note that, for a rigid solid, the deformation is absent, thusũ x =ũ y = 0, which removes the Dv (i) xũ y terms in the tangential velocity continuity equations (2.16) and (2.20), eliminating the instability. Interestingly, from table 3, the decay rate of the second (stable) eigenvalue is less affected by the absence of the terms Dv (i) xũ y , compared to the unstable eigenvalue, while the decay rates of the remaining two eigenvalues remain unaffected even after removal of the interaction terms (Dv (i) xũ y ). This implies that the unstable mode is destabilised due to the energy exchange at the fluid-DS interface while other modes remain unaffected or are very weakly affected.
The role of μ r in the unconditionally unstable region is shown in figure 7(a). In contrast to the case of high a, for sufficiently low a the coupling established between the fluids through the DS is strong enough to make the coupled system unstable at non-zero Γ . However, the unconditionally unstable region comes at the 'price' of a low growth rate at low Γ , which will be further discussed in the subsequent section.
The important role of fluid 2 in the induced instabilities can be illustrated with the help of figure 7(c). For b = 0.1, the Γ c curve is very similar to the curve for a Couette flow past an elastic solid, since Γ c decreases with increasing a. This implies that a minimum a μ r = 0.5 μ r = 0.75 μ r = 0.9 μ r = 1 μ r = 5 a V r = -1 V r = 0 V r = 0.5 10 -1 10 0 10 1 10 1 10 1 10 0 10 -1 10 -1 10 -2 10 -3 10 0 a Γ thickness of fluid 2 is necessary for the destabilising coupling between the two fluids via the DS. For b = 1, this requirement is fulfilled and the coupling has a strong effect on the viscous mode instability. Furthermore, increasing b leads to an increase in the minimum value of a, at which point the system becomes unconditionally unstable. Figure 7(b) illustrates the effect of the variation in V r , which represents the shear force acting on fluid 2, compared to fluid 1. The figure shows that the co-flow configuration is more unstable than counter-flow. As explained in the preceding discussion, the terms Dv (i) xũ y in the tangential velocity continuity equations (2.16) and (2.20) are the main cause of the instability predicted in the present work. The influence of the coupling terms increases with an increase in the DS deformation, which is, in turn, induced by the shear force exerted by the flow. In the counter-flow configuration, the shear force acting on the DS due to both fluids is in the opposite direction, reducing the deformation in the DS and, hence, stabilising the disturbances. Conversely, under co-flow, the shear forces exerted by the fluids reinforce each other, increasing the coupling and resulting in de-stabilisation. The effect of varying the viscosity ratio, μ r = μ 2 /μ 1 on the critical parameters is shown in figure 7(a). For μ r = 1, the unconditionally unstable region does not exist, suggesting a strongly stabilising effect imparted by the viscosity of fluid 2. Through computations of similar Γ c curves for various values of μ r in the range [0.5, 1], the critical value of μ r was found to be ∼ 0.98 for the existence of the unconditionally unstable region. Furthermore, the curve corresponding to μ r = 0.75 illustrates that an increase in viscosity of fluid 2 can push the unconditionally unstable region to lower values of a. Figures 7(a), 7(b) and 7(c) also show that in certain parametric regimes, Γ c ∼ O(1) for which the linear elastic model is not applicable (Gaurav Shankar 2009, 2010Patne & Shankar 2018). Instead, a material frame-invariant model such as the neo-Hookean or Mooney-Rivlin model is applicable. However, in applications such a high value of Γ is not realisable owing to the practical and application restrictions. The results for such a high Γ have been presented here only for the sake of completion.
In figures 6, 7(a), 7(b) and 7(c), the system is unstable for any non-zero value of Γ and k c → 0. However, for such an unconditionally unstable region, the growth rate of the disturbances also decreases with decreasing k and Γ , as illustrated in figure 7(d).
Owing to the small growth rate, this regime of the instability may not actually be observable experimentally. We note that the low a regime may be relevant in processes involving thin elastic films such as membranes or lipid bilayers. Furthermore, to induce mixing in a practical setting, the growth rate should be high, such that it may lead to the rapid movement of fluid particles (Verma & Kumaran 2012Bodiguel et al. 2015). Figure 7(d) shows that, for k = 0.005 and Γ = 10 −10 , the characteristic growth time of the disturbances is ∼10 7 s, while for k = 0.5 and Γ = 10 −4 , it is ∼10 2 s, which is realistically achievable experimentally and in practice. In applications, owing to the practical constraints, the length of the DS will be restricted. Accordingly, for a DS of length L, the smallest disturbance wavenumber would be 2π/L. An estimate to observe the coupling instability predicted here in the experiments or practical applications can be obtained as follows. Consider the membrane separation process for which the dimensionless parameters from table 1 are a ∼ O(10 −3 − 10 −6 ), b ∼ O(10 −3 − 10 −6 ), V r ∼ O(0 − 10) and μ r ∼ O(0.1 − 10). To detect the coupling mode instability, assume μ r ∼ 0.5 and V r = 0 so that the figure 7(d) can be readily used. Then the instability exists for k < 0.9 and arbitrary Γ values. For membrane separation processes, length L ∼ O(10 −2 − 10) m. Thus, depending on L, the minimum k will be selected by the system. For example, if L ∼ O(1), then from figure 7(d), the temporal coupling instability will exhibit a growth rate of O(10 −2 ) s −1 which is possible to be observed in the experiments. Also, the dimensionless speed of the plate needed to destabilise the disturbances having k ∼ 0.8 is Γ ∼ 10 −4 which corresponds to V 1 ∼ 10 −2 − 10 −4 m s −1 , for a membrane of shear modulus G ∼ 10 5 Pa and channel height R ∼ 10 −4 − 10 −6 m which again falls within a practical parametric regime. Please note that, if the interfacial tension is sufficiently small, then the coupling instability predicted here can exist for k > 1, in which case the minimum L required to observe instability can become as low as ∼ 0.1.
Further analysis at k → 0 reveals that, for any given combination of the parameters, k c / = 0 since for k → 0, only one non-trivial eigenvalue exists, which is unconditionally stable. To conclude, for the low a regime, despite exhibiting instability at small k and Γ , due to practical restrictions and small growth rate, such instability may not be of physical relevance. However, the same instability might play a significant role in inducing mixing at sufficiently high k and Γ .

Instability characteristics at finite Re
The results discussed in § 3.1.1 are for a vanishing Reynolds number, Re. However, physical restrictions dictate that Re may be small but of a finite value. Furthermore, the previous study of Thaokar & Kumaran (2002) only considered the creeping-flow limit in their analysis, neglecting the influence of inertia on the instability. Motivated by these two reasons, the effect of the variation in Re on the instabilities is studied in the present section.
In the creeping-flow limit, as shown in table 2, only four eigenvalues exist for arbitrary values of the parameters. However, for non-zero Re there is an unlimited number of eigenvalues, dependent on the number of collocation points used to discretise the differential operator. For N = 40 collocation points, the eigenspectrum is shown in figure 8. The eigenvalues, other than those present in the creeping-flow limit, are the shear waves generated at the fluid-DS interfaces. Interestingly, these shear waves are stable for Re < 1. For Re > 1, these waves may give rise to multiple unstable modes (Gaurav Shankar 2009, 2010Patne & Shankar 2018). In the present system, as discussed in § 3.1.1, Re < 1, thus shear waves are found to be stable. This implies that, unlike in the viscous mode, the shear waves remain unaffected by the presence of fluid 2 at low Re. As a consequence of the stability of the shear waves for Re < 1, the viscous mode determines the stability of the system even in the presence of finite inertia. Figure 9 illustrates the effect of the increase in Re on the critical parameters. Increasing inertia has a slightly stabilising effect on the viscous mode. Thus, the results presented for Re = 0 in § 3.1.1 are generally applicable for Re < 1, with a minor quantitative adjustment.

Spatio-temporal stability analysis (k i
When performing a temporal stability analysis (see § 3.1), the frequency of the disturbances, ω, is taken as a complex number and the wavenumber, k, is a real number (Drazin & Reid 1981). Conversely, for a spatial stability analysis, used to study the evolution of disturbances in space, ω is taken as a real number and k as a complex number (Drazin 2002). In a spatio-temporal stability analysis, both ω and k are treated as complex numbers (Huerre & Monkewitz 1990;Schmid & Henningson 2001) and instability is classified as either an absolute or convective instability. This section is aimed at understanding the effect of the fluid-solid coupling on the spatio-temporal instabilities, in the creeping-flow limit. The classification of flows as absolutely or convectively unstable was first proposed by Briggs (1964) in the context of plasma physics. He further developed the methodology to investigate the spatio-temporal instabilities by progressive moving of the Laplacian contour in the complex frequency plane and Fourier contour in the wavenumber plane. According to Briggs (1964), absolute instability signifies the growth of disturbances in both upstream and downstream directions, while convective instability implies that the disturbances develop in the downstream direction from the source of the disturbances. Thus, convective instability will decay at any fixed position in space if provided sufficient time (Huerre & Monkewitz 1990), resulting in mixing only in the downstream direction, while absolute instability will induce mixing both up-and downstream. This can be illustrated by considering the response of a given base velocity profile to an impulse excitation at asymptotically long times (Huerre & Monkewitz 1990), where the obtained response is used to determine whether the flow is absolutely or convectively unstable.
The absolute instability is absent in a plane Couette flow past a rigid surface since the flow is temporally stable (Drazin 2002) -for the existence of an absolute or convective instability, the flow must be temporally unstable (Huerre & Monkewitz 1990;Schmid & Henningson 2001). However, the plane Couette flow past a plate lined by a DS exhibits an absolute instability since the fluid-DS coupling gives rise to the temporal instability (Patne & Shankar 2017). The present flow geometry also possesses a deformable boundary but it differs from the plane Couette flow past a DS due to the presence of fluid 2, which may lead to modification of the results predicted by Patne & Shankar (2017). The assumption of creeping flow in the present section is due to the considerable simplification in the numerical determination of the saddle and cusp points. Also, as shown in § 3.1.2, Re has a negligible effect on the instabilities for Re < 1.
To illustrate the absolute and convective instabilities, consider a dispersion relation of the form ω = f (k), where f (k) is a continuous and differentiable function of k. For the absolute instability to exist, the group velocity of the disturbances must vanish for at least one value of k, so that ∂ω/∂k = 0 (Huerre & Monkewitz 1990;Schmid & Henningson 2001). However, the determined saddle point must also obey the causality principle for the existence of the absolute instability (Huerre & Monkewitz 1990). To obtain the sufficient condition, equation ∂ω/∂k = 0 is solved for k, which gives the saddle points of the dispersion relation. If the saddle point is of first order in the k-plane (denoted by k 0 ) and ω 0 is the value of ω at the saddle point k 0 , then a local Taylor expansion about this point gives (ω − ω 0 ) ∼ (k − k 0 ) 2 . This shows that the mapping from the k-plane to the ω-plane is characterised by angle doubling (i.e. phase doubling) provided that the saddle point is of first order. Here, the k-and ω-planes refer to the k r − k i and ω r − ω i planes, respectively.
For realistic dispersion relations, the solution of the equation ∂ω/∂k = 0 gives several roots and thus saddle points. To determine the genuine saddle point, obeying the causality principle (which stipulates that the cause does not precede the effect), the method of Briggs (1964) is needed. However, finding the saddle point in the k-plane and cusp point in ω-plane using the method of Briggs (1964) becomes a cumbersome mathematical and numerical exercise for realistic dispersion relations. A simpler alternative is the method of Kupfer, Bers & Ram (1987), in which for prediction of an absolute instability only the formation of the cusp point in the ω-plane is necessary.
Use of the method presented by Kupfer et al. (1987) is illustrated in figure 10, where the cusp point forms in the ω-plane at ω 0 = 0.1142 + 0.01734i. To obtain the cusp point, the temporal stability curve is first obtained in the ω-plane by varying k r and setting k i = 0. Next, k i is given a negative value so that the disturbances can become spatially unstable. Slowly lowering k i eventually forms a cusp point at k i = −0.24125, as shown in figure 10. To verify the genuine character of the cusp point, a straight line (not shown in figure 10) is drawn from the cusp point parallel to the ω i axis. The intersections of the drawn straight line and the temporal stability curve are then counted. For example, in figure 10 the count is one -an odd number -thus the formed cusp point is genuine. It must be noted that if the count is even then it is an evanescent cusp point. The existence of a genuine cusp point implies the presence of the absolute instability, while that of an evanescent mode signifies a spurious cusp point (Yeo, Khoo & Zhao 1996, 1999. According to the method of Kupfer et al. (1987), a cusp point in the ω-plane must correspond to a saddle point in the k-plane. To obtain the saddle point in the k-plane, a matrix of ω r and ω i values is generated by varying both k r and k i and (using MATLAB built-in functions) isocontours are obtained as shown in figure 11. The cusp point (ω 0 ) of figure 10 and saddle point (k 0 ) of figure 11 represent a local mapping of the type ω − ω 0 ∼ (k − k 0 ) 2 , which shows that the saddle point in the k-plane is of first order.
Computation of the cusp point is a time-consuming procedure since an analytical form of the dispersion relation is not available. Thus, unlike for the temporal stability analysis, critical parameters are not determined. Instead, parameters are varied in figures 12 and 13(a) with respect to those used to obtain figure 10  is noted and then used in determining the stabilising or destabilising role of the parameters. In figure 12(a), Γ is decreased with respect to figure 10 and leads to the absence of the absolute instability, implying a stabilising effect of decreasing Γ . This also shows that for the existence of the absolute instability, similar to the temporal instability, there exists a critical value of Γ at which the flow undergoes a transition from a convectively unstable to an absolutely unstable state. Figure 12  thickness of the DS, where increasing a by two orders of magnitude compared to figure 10 leads to the disappearance of the absolute instability. This also shows that even when the temporal analysis exhibited an unconditionally unstable region at low a, it is absent for the absolute instability.
In figure 12(c), compared to figure 10, b is decreased from b = 2 to b = 1, resulting in a decreased growth rate of the absolute instability (ω i0 ) by one order of magnitude, demonstrating the strong stabilising effect of decreasing b. The effect of varying μ r on the instability is illustrated in figure 12(d), where μ r is increased from μ r = 0.1 (figure 10) to μ r = 0.5. As seen with the temporal instability, increasing μ r has a stabilising effect on the absolute instability.
Lastly, figure 13(a) shows that the counter-flow configuration is favourable for inducing the absolute instability, illustrated by the fact that ω i0 is higher for V r = −1 than V r = 0 (figure 10). The absolutely unstable flow implies that the disturbances must spatio-temporally grow both in the downstream and upstream directions. We suspect that the counter-flow configuration is more easily amenable to the spatio-temporal growth of the disturbances in the upstream direction since fluid 2 is exerting a shear stress on the DS in the upstream direction due to the flow in negative x-direction. Such a shear stress can then help amplify the disturbances in the negative x-direction. Furthermore, fluid 2 can also convect the disturbances in the negative x-direction under the counterflow configuration, which is essential for the spatial convection of the disturbances. These could be the reasons that the counterflow is more prone to an absolute instability.
The preceding results are concerned with the existence of the absolute instability at low a. However, the absolute instability also exists at higher a, provided that Γ is also sufficiently high, as shown in figure 13(b). To conclude, by carefully adjusting the properties of the fluids and the DS, it is possible to introduce absolute instability in the system, which can in turn be used potentially to enhance mixing.
The parametric regime to observe the absolute instability in the experiments can be obtained as follows. From figure 10, the absolute instability exists for the parameters a = 10 −4 , b = 2, μ r = 0.1, V r = 0, Γ = 10 −3 , Σ 1 = 1 and Σ 2 = 1 which corresponds to V 1 ∼ 10 −1 − 10 −3 m s −1 , for a membrane of shear modulus G ∼ 10 5 Pa and channel height R ∼ 10 −4 − 10 −6 m which is achievable in the practical applications. Furthermore, from figure 10 and 11, the spatial and temporal growth rates of the absolute instability at the above parameters are ω i = 0.01734 and k i = 0.24125, which seem detectable, experimentally. Using k r value from figure 11, the minimum length of the DS necessary to observe the above instability is of O(1) m.

Conclusions
The present study examined the linear temporal and spatio-temporal stability of simple shear flows coupled by an elastic, deformable solid layer. The fluids are Newtonian while the DS is described using linear elasticity. Such a fluid-DS-fluid coupling can be encountered in microfluidics, membrane separation and formation as well as in physiological systems such as cells and the heart. The present study illustrates the effect of the coupling between the fluids, as mediated by a finite-thickness DS, on the elastohydrodynamic instability and the existence of absolute instability.
At a sufficiently high dimensionless thickness of the DS (a), a viscous instability is predicted -in agreement with the plane Couette flow past a plate lined by a DS. Under these conditions, the two flows on either side of the DS are hydrodynamically decoupled. However, at low a, the flows interact via the DS, which affects the viscous mode of the instability. The predicted instability arises from the energy exchange between the fluids and the DS via the tangential velocities at the interface. The interaction lowers the value of Γ c (Γ represents the shearing motion) required to destabilise the finite wave disturbances, eventually leading to an unconditionally unstable region. The existence of the unconditionally unstable region strongly depends on the thickness and viscosity of fluid 2, i.e. in dimensionless terms b and μ r , respectively. While this unstable region is present at low Γ and k (the wavenumber), it is accompanied by a very small growth rate of the perturbations. Such a low value of Γ and practical restrictions on k define a parametric limit within which the predicted instabilities may be observed in experiments and applications.
To further understand the role of the DS in triggering and controlling potential mixing in practical applications, a spatio-temporal stability analysis was carried out. The analysis revealed the existence of the absolute instability in the present system, however, at much higher Γ compared to that required for triggering the temporal instability. At low a, since the growth rate of the disturbances is small, the existence of the absolute instability can be highly useful since it will enhance mixing. Finally, in practical applications, the DS can be of a length comparable to the widths of the fluids. In such cases, a full global stability analysis is required and may give rise to new types of instabilities, presenting a possible next step in understanding the effect of the fluid-fluid coupling via a deformable solid on elastohydrodynamic instabilities.