1 Introduction
Two-fluid flows display interesting instabilities due to viscosity and density contrasts between the fluids. Differences in these properties across the flow often exist simultaneously, but our interest is in isolating the effects of viscosity contrasts alone. The reverse case, of contrasting density but constant viscosity, has been far more widely studied in geophysical and other contexts. Instabilities due to viscosity variation too have been investigated by several authors (see e.g. Yih Reference Yih1967; Preziosi, Chen & Joseph Reference Preziosi, Chen and Joseph1989; Chen & Meiburg Reference Chen and Meiburg1996; Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Joseph et al. Reference Joseph, Bai, Chen and Renardy1997; Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999; Balasubramaniam et al. Reference Balasubramaniam, Rashidnia, Maxworthy and Kuang2005; d’Olce et al. Reference d’Olce, Martin, Rakotomalala, Salin and Talon2008; Talon, Goyal & Meiburg Reference Talon, Goyal and Meiburg2013; John et al. Reference John, Oliveira, Heussler and Meiburg2013; Naraigh et al. Reference Naraigh, Valluri, Scott, Bethune and Spelt2014) for both immiscible and miscible fluids. An extensive discussion of instability associated with such flows can be found in a recent review by Govindarajan & Sahu (Reference Govindarajan and Sahu2014).
By conducting a linear-stability analysis, Yih (Reference Yih1967) was the first to demonstrate that immiscible fluid layers (with a sharp interface between the two fluids) in shear flow are unstable to infinitesimally small long-wave disturbances at any Reynolds number. Since then, instability in the context of sharp interfaces has been investigated by many researchers (e.g., Hooper & Boyd Reference Hooper and Boyd1983; Hooper Reference Hooper1985; Valluri et al. Reference Valluri, Naraigh, Ding and Spelt2010), and short-wave instabilities were found as well. The mechanism of this short-wave instability was provided by Hinch (Reference Hinch1984).
Miscible flows are different from immiscible flows in an important way. The diffusivity, characterised by the inverse of the Schmidt number, is among the factors that play an important role. Govindarajan (Reference Govindarajan2004) investigated three-layer Poiseuille channel flows (wherein two miscible fluids are separated by a mixed region) and showed that at higher Schmidt numbers these flows go unstable at lower Reynolds numbers. She found that when the more (less) viscous fluid occupies the near-wall regions, the flow is significantly destabilised (stabilised). These effects are accentuated by an increase in viscosity contrast. For pipe flow when the viscosity ratio is large, Selvam et al. (Reference Selvam, Merk, Govindarajan and Meiburg2007) found that the flow can be destabilised even in the opposite scenario, i.e. when the less viscous fluid is near the wall. Ern, Charru & Luchini (Reference Ern, Charru and Luchini2003) studied the influence of diffusion and mixed layer thickness in a miscible two-fluid Couette flow at a Reynolds number of less than one, and showed that diffusivity has a non-monotonic effect on the growth rate of the disturbance. They reported that flows at intermediate Schmidt numbers can be more unstable than flows at either very low or very high Schmidt numbers. They also found regimes at low Reynolds numbers where miscible flows are more unstable than those interfacial flows (with no diffusion across the sharp interface). The geometry and flow considered in the present work is the same as that of Talon & Meiburg (Reference Talon and Meiburg2011). They found that in the limit
$\mathit{Re}\rightarrow 0$
, instability can be triggered by four different types of modes depending on the interface location. A mechanism of these instabilities in the Stokes flow regime was also provided based on perturbation of the concentration relative to the interface. The authors distinguish this instability from the inertial mechanisms of both Hinch (Reference Hinch1984) and Govindarajan (Reference Govindarajan2004). Talon & Meiburg (Reference Talon and Meiburg2011) also remarked that the instability observed in their study is similar to that of Ern et al. (Reference Ern, Charru and Luchini2003) in the Couette flow. Instability at high Reynolds number of this miscible flow and its dynamics in the nonlinear regime has not been investigated yet. In the present study, we also investigate the difference between interfacial instabilities (without surface tension) with a viscosity jump across the interface, and instabilities at high Schmidt number (poor diffusivity) and at high Reynolds number in pressure-driven two-layer miscible channel flow, by performing linear stability analyses and direct numerical simulations.
Important in our discussion will be the location of the critical layer (the layer at which the phase speed of the disturbance is close to the mean streamwise velocity, and a major portion of the kinetic energy production takes place). In miscible flow, when this layer overlaps the viscosity-stratified layer, the dominant balance at the lowest order is changed (Govindarajan Reference Govindarajan2004), so an additional ‘overlap’ mode of instability can occur. When the two layers are well separated, the Tollmien–Schlichting mode of instability is the most likely. In immiscible flows too, both scenarios can occur, i.e. the interface may or may not coincide with the critical location of the dominant disturbance.
The rest of the paper is organised as follows: the mathematical formulation of the linear stability equations for miscible and immiscible flows are presented in §§ 2 and 3, respectively. The budget for disturbance kinetic energy is formulated in § 4. The results are discussed in § 5, and concluding remarks are given in § 6.
2 Formulation: two miscible fluids with a mixed layer in between
The linear stability analysis and direct numerical simulation of a two-layer channel flow made up of two miscible, Newtonian and incompressible fluids of equal density and different viscosities is considered. The Cartesian coordinate system
$(x,y,z)$
is used to formulate the problem, where
$x$
,
$y$
and
$z$
denote the coordinates in the horizontal, the vertical and the spanwise directions, respectively. As shown in figure 1(a), the top (fluid ‘1’ of dynamic viscosity
${\it\mu}_{1}$
) and bottom (fluid ‘2’ of dynamic viscosity
${\it\mu}_{2}$
) fluids occupy the regions
$h+q/2\leqslant y\leqslant H$
and
$0\leqslant y\leqslant h-q/2$
, respectively, where
$q$
is the mixed layer thickness. The channel walls are located at
$y=0$
and
$y=H$
, and periodic boundary conditions are imposed in the spanwise direction. The viscosity variation occurs due to the spatially varying magnitude of a scalar (
$s$
), which could be, for example, the concentration of a solute, or temperature. Without loss of generality, the base-state concentration
$s_{0}$
is taken to be 1 in the bottom layer, and 0 in the top layer, and it varies from 1 to 0 in the mixed layer. Thus,
$q=H$
and
$0$
represent complete stratification and a sharp interface, respectively. The latter is shown in figure 1(b).
Figure 1. Schematic of the flows. (a) Miscible case: fluids ‘2’ and ‘1’ occupy the bottom (
$0\leqslant y\leqslant h-q/2$
) and top (
$h+q/2\leqslant y\leqslant H$
) layers, respectively. The two fluids are separated by a mixed layer of uniform thickness
$q$
. (b) The corresponding immiscible case.
The viscosity,
${\it\mu}$
, is modelled as an exponential function of the scalar
$s$
:
$$\begin{eqnarray}{\it\mu}={\it\mu}_{1}\text{exp}\left(sR_{s}\right),\end{eqnarray}$$
where
$R_{s}\left(\equiv \mathit{ln}\left({\it\mu}_{2}/{\it\mu}_{1}\right)\right)$
is the log-mobility ratio of the scalar. The following scaling is employed to render the governing equations dimensionless:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (x,y,z,q,h)=H(\widetilde{x},\widetilde{y},\widetilde{z},\widetilde{q},\widetilde{h}),\quad t=\frac{H^{2}}{Q}\widetilde{t},\quad (u,v)=\frac{Q}{H}(\widetilde{u},\widetilde{v}),\\ \displaystyle p=\frac{{\it\rho}Q^{2}}{H^{2}}\widetilde{p},\quad {\it\mu}=\widetilde{{\it\mu}}{\it\mu}_{1},\end{array}\right\}\end{eqnarray}$$
where the tildes designate dimensionless quantities:
$Q$
denotes the total volume flow rate per unit distance in the spanwise direction;
$u$
,
$v$
and
$w$
are the velocity components in the
$x$
,
$y$
and
$z$
directions, respectively;
$p$
denotes pressure;
${\it\rho}$
is the constant density and
$t$
is time. The dimensionless governing equations (after dropping the tildes) are given by
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$
$$\begin{eqnarray}\left[\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}\right]=-\boldsymbol{{\rm\nabla}}p+\frac{1}{\mathit{Re}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}})\right],\end{eqnarray}$$
$$\begin{eqnarray}\frac{\partial s}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s=\frac{1}{\mathit{Sc}\mathit{Re}}{\rm\nabla}^{2}s,\end{eqnarray}$$
where
$\boldsymbol{u}$
is the velocity vector, and
$\mathit{Re}\,(\equiv {\it\rho}Q/{\it\mu}_{1})$
and
$\mathit{Sc}\,(\equiv {\it\mu}_{1}/{\it\rho}\mathscr{D})$
are the Reynolds number and Schmidt number, respectively, wherein
${\mathcal{D}}$
is the diffusion coefficient of the scalar.
2.1 Base state
The base state corresponds to a steady, parallel, fully developed flow, i.e.
$U=U(y),V=W=0$
and
$P$
is linear in
$x$
. Here, the base-state quantities are designated by upper-case letters for the flow variables, and by the subscript
$0$
for viscosity and
$s$
. In order to make the concentration of the scalar continuous up to the second derivative at
$y=h-q/2$
and
$y=h+q/2$
, the mean scalar
$s_{0}(y)$
is chosen to be fifth-order polynomial in the mixed layer (Malik & Hooper Reference Malik and Hooper2005):
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle s_{0}=1,\quad 0\leqslant y\leqslant h-q/2,\\ \displaystyle s_{0}=\mathop{\sum }_{i=1}^{6}a_{i}y^{i-1},\quad h-q/2\leqslant y\leqslant h+q/2,\\ \displaystyle s_{0}=0,\quad h+q/2\leqslant y\leqslant 1,\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where the
$a_{i}~(i=1,6)$
are given by
$$\begin{eqnarray}\displaystyle \hspace{-19.20007pt}\left.\begin{array}{@{}c@{}}\displaystyle a_{1}=\frac{(h+q/2)^{3}}{q^{5}}[6(h-q/2)^{2}-3(h-q/2)q+q^{2}],\quad a_{2}=-\frac{30(h-q/2)^{2}(h+q/2)^{2}}{q^{5}},\\ \displaystyle a_{3}=\frac{60(h-q/2)}{q^{5}}(h+q/2)h,\quad a_{4}=-\frac{10}{q^{5}}[6(h-q/2)^{2}+6(h-q/2)q+q^{2}],\\ \displaystyle a_{5}=\frac{30}{q^{5}}h\quad \text{and}\quad a_{6}=-\frac{6}{q^{5}}.\end{array}\right\}\hspace{-20.00003pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
We have confirmed that results indistinguishable from the present ones are obtained by using any other sufficiently smooth profile, such as the error function or the hyperbolic tangent. For the parameter range of the present study (high Reynolds numbers and high Péclet numbers) the mixed layer diffuses very slowly, with a divergence angle of the order of
$Pe^{-1}$
. Thus the assumption of locally parallel flow and the ensuing use of a constant thickness mixed layer are extremely reasonable in this context, with errors of
$O[Pe^{-1}]$
. A brief description on the validity of the parallel flow assumption is provided in appendix A. The assumption will also be justified later in the form of comparisons with direct numerical simulations, where no such assumption is made.
The base-state streamwise velocity profile
$U(y)$
is obtained by solving the steady, fully developed version of (2.4) using no-slip and no-flux conditions at the wall and the centreline of the channel, respectively, i.e.,
$$\begin{eqnarray}\mathit{Re}\left(\frac{\text{d}P}{\text{d}x}\right)=\left({\it\mu}_{0}U^{\prime }\right)^{\prime },\end{eqnarray}$$
where
${\it\mu}_{0}=\text{e}^{\left(R_{s}s_{0}\right)}$
and the prime represents differentiation with respect to
$y$
. The non-dimensional pressure gradient
$\text{d}P/\text{d}x$
is fixed by using
$\int _{0}^{1}U\,\text{d}y=1$
.
Figure 2. Typical base steady-state profiles of (a)
${\it\mu}_{0}$
and (b)
$U^{\prime \prime }$
for different values of
$R_{s}$
. The other parameters are chosen as
$h=0.25$
and
$q=0.1$
.
In figure 2(a,b) we show, respectively, typical profiles of the base-state viscosity and of the second derivative of the mean velocity for different values of
$R_{s}$
. We choose to show the second derivative rather than the velocity profile itself because it demonstrates that the velocity profile changes slope very rapidly in the mixed layer. Moreover, the case of
$R_{s}=1$
is seen to contain a point of inflexion, which indicates a tendency for inviscid instability.
2.2 Linear stability analysis
The temporal linear stability of the base flow given by (2.6)–(2.8) using a normal mode analysis considering two-dimensional perturbations is investigated. For a single fluid flow, Squire’s theorem (Squire Reference Squire1933) states that every unstable three-dimensional disturbance is associated with an equally unstable two-dimensional disturbance at a lower value of Reynolds number. We assume that in our stratified flow too, two-dimensional disturbances go unstable at lower Reynolds numbers than three-dimensional ones. We therefore study the linear instability to two-dimensional perturbations. Our direct numerical simulations confirm that in the range we study, two-dimensional disturbances are the first to go unstable. The flow variables are split into base-state quantities and two-dimensional perturbations (designated by a hat):
$$\begin{eqnarray}(u,v,p,s)(x,y,t)=(U(y),0,P,s_{0}(y))+(\hat{u} ,\hat{v},\hat{p},{\hat{s}})(y)\text{e}^{\text{i}({\it\alpha}x-{\it\omega}t)},\end{eqnarray}$$
and the perturbation viscosity is given by
$$\begin{eqnarray}\hat{{\it\mu}}=\frac{\text{d}{\it\mu}_{0}}{\text{d}s_{0}}{\hat{s}},\end{eqnarray}$$
where
$\text{i}\equiv \sqrt{-1}$
,
${\it\alpha}$
and
${\it\omega}\,(\equiv {\it\alpha}c)$
are the wavenumber and frequency of the disturbance, respectively, and wherein
$c$
is the phase speed of the disturbance. In temporal stability analysis,
${\it\alpha}$
and
${\it\omega}$
are treated as real and complex quantities, respectively, whereas both are complex in spatio-temporal analysis (e.g. see Schmid & Henningson Reference Schmid and Henningson2001). We conduct the former here, where a given mode is unstable if
${\it\omega}_{i}>0$
, stable if
${\it\omega}_{i}<0$
and neutrally stable if
${\it\omega}_{i}=0$
;
${\it\omega}_{i}$
being the imaginary part of
${\it\omega}$
.
Following a standard approach by substituting (2.9) into (2.3)–(2.5), subtracting the base-state equations and subsequently linearising and eliminating the pressure perturbation, we obtain the following linear stability equations (Govindarajan Reference Govindarajan2004), with the hat notation suppressed:
$$\begin{eqnarray}\displaystyle & & \displaystyle \text{i}{\it\alpha}\mathit{Re}[({\it\psi}^{\prime \prime }-{\it\alpha}^{2}{\it\psi})(U-c)-U^{\prime \prime }{\it\psi}]={\it\mu}_{0}({\it\psi}^{iv}-2{\it\alpha}^{2}{\it\psi}^{\prime \prime }+{\it\alpha}^{4}{\it\psi})\nonumber\\ \displaystyle & & \displaystyle \quad +2{\it\mu}_{0}^{\prime }({\it\psi}^{\prime \prime \prime }-{\it\alpha}^{2}{\it\psi}^{\prime })+{\it\mu}_{0}^{\prime \prime }({\it\psi}^{\prime \prime }+{\it\alpha}^{2}{\it\psi})+U^{\prime }({\it\mu}^{\prime \prime }+{\it\alpha}^{2}{\it\mu})+2U^{\prime \prime }{\it\mu}^{\prime }+U^{\prime \prime \prime }{\it\mu},\qquad\end{eqnarray}$$
$$\begin{eqnarray}\text{i}{\it\alpha}\mathit{Sc}\mathit{Re}[(U-c)s-{\it\psi}{s_{0}}^{\prime }]=(s^{\prime \prime }-{\it\alpha}^{2}s),\end{eqnarray}$$
wherein the amplitude of the velocity disturbances is re-expressed in terms of a streamfunction [
$(\hat{u} ,\hat{v})=({\it\psi}^{\prime },-\text{i}{\it\alpha}{\it\psi})$
].
Solutions of these equations are obtained subject to the following boundary conditions at both walls
$$\begin{eqnarray}\displaystyle {\it\psi}={\it\psi}^{\prime }=s^{\prime }=0\quad \text{at }y\pm 1. & & \displaystyle\end{eqnarray}$$
Equations (2.11)–(2.12) along with the boundary conditions (2.13) constitute an eigenvalue problem, which is solved using the public domain software LAPACK. A Chebyshev spectral collocation is used to discretise the domain. Due to the presence of large gradients in the viscosity-stratified region, a large number of grid points are required in this region. For this we use the stretching function proposed by Govindarajan (Reference Govindarajan2004):
$$\begin{eqnarray}y_{j}=\frac{a}{\sinh (by_{0})}[\sinh \left\{(y_{c}-y_{0})b\right\}+\sinh (by_{0})],\end{eqnarray}$$
where
$y_{j}$
are the locations of the grid points,
$a$
is the mid-point of the stratified layer,
$y_{c}$
is a Chebyshev collocation point,
$$\begin{eqnarray}y_{0}=\frac{0.5}{b}\mathit{ln}\left[\frac{1+(e^{b}-1)a}{1+(\text{e}^{-b}-1)a}\right],\end{eqnarray}$$
and
$b$
is the degree of clustering;
$b=8$
is taken in this present study. The above formulation gives an accuracy of at least five decimal places in the range of parameters used.
3 Formulation: two immiscible fluids separated by a sharp interface
3.1 Base state
For pressure-driven flow of two immiscible fluids separated by a sharp interface (shown in figure 1 b), the base-state velocity profile is given by
$$\begin{eqnarray}\displaystyle & \displaystyle U_{1}=\frac{1}{2}\frac{\left(\text{d}P/\text{d}x\right)^{-1}}{\text{e}^{R_{s}}}\left[\frac{\text{d}P}{\text{d}x}y+c_{3}\right]^{2}+c_{4}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle U_{2}=\frac{\text{d}P}{\text{d}x}\frac{y^{2}}{2}+c_{1}y+c_{2}, & \displaystyle\end{eqnarray}$$
with subscripts
$1$
and
$2$
denoting the upper and lower layers, respectively. We obtained (3.1) and (3.2) by integrating the steady, fully developed dimensionless Navier–Stokes equations. Taking the undisturbed height of the interface to be
$h_{0}$
, the pressure gradient,
$\text{d}P/\text{d}x$
, and the integration constants,
$c_{1}$
,
$c_{2}$
,
$c_{3}$
and
$c_{4}$
, are obtained by solving the following simultaneous equations, which correspond to no-slip conditions at the walls and balance of the tangential component of the stress at the interface.
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\left(\text{d}P/\text{d}x\right)^{-1}}{2\text{e}^{R_{s}}}\left\{\left[\frac{\text{d}P}{\text{d}x}h_{0}+c_{3}\right]^{2}-c_{3}^{2}\right\}-\frac{1}{2}\frac{\text{d}P}{\text{d}x}\left({h_{0}}^{2}-1\right)-c_{1}\left(h_{0}-1\right)=0,\\ \displaystyle c_{3}=c_{1},\quad c_{2}=-\frac{1}{2}\frac{\text{d}P}{\text{d}x}-c_{1},\quad c_{4}=-\frac{\left(\text{d}P/\text{d}x\right)^{-1}}{2\text{e}^{R_{s}}c_{1}^{2}}.\end{array}\right\}\end{eqnarray}$$
The pressure gradient,
$\text{d}P/\text{d}x$
, is obtained from the constant volumetric flow rate condition, i.e.,
$$\begin{eqnarray}\int _{0}^{h_{0}}U_{1}\,\text{d}y+\int _{h_{0}}^{1}U_{2}\,\text{d}y=1.\end{eqnarray}$$
3.2 Linear stability analysis
We also examine the linear stability of the base state, obtained by solving (3.1) and (3.2), to infinitesimal, two-dimensional disturbances. Each flow variable is expressed as the sum of a base state and a two-dimensional perturbation,
$$\begin{eqnarray}(\tilde{u} _{i},\tilde{v}_{i},\tilde{P}_{i})(x,y,t)=\left[U_{i}(y),0,P_{i}\right]+\left(\hat{u} _{i},\hat{v}_{i},\hat{p}_{i}\right)(x,y,t),\end{eqnarray}$$
with
$i=1,2$
. Similarly the height
$h$
of the interface can be expressed as
$$\begin{eqnarray}h(x,y,t)=h_{0}+{\hat{h}}\quad \text{with }{\hat{h}}(x,t)=\tilde{h}\text{e}^{\text{i}({\it\alpha}x-{\it\omega}t)}.\end{eqnarray}$$
Substituting (3.5) and (3.6) into the governing equations and following the same procedure as before yields the following linear stability equations. In the lower layer:
$$\begin{eqnarray}\text{i}{\it\alpha}\mathit{Re}[\{{v_{1}}^{\prime \prime }-{\it\alpha}^{2}v_{1}\}(U_{1}-c)-U_{1}^{\prime \prime }v_{1}]=v_{1}^{\prime \prime \prime \prime }-2{\it\alpha}^{2}v_{1}^{\prime \prime }+{\it\alpha}^{4}v_{1}.\end{eqnarray}$$
In the upper layer:
$$\begin{eqnarray}\text{i}{\it\alpha}\mathit{Re}[\{{v_{2}}^{\prime \prime }-{\it\alpha}^{2}v_{2}\}(U_{2}-c)-U_{2}^{\prime \prime }v_{2}]=\text{e}^{R_{s}}[v_{2}^{\prime \prime \prime \prime }-2{\it\alpha}^{2}v_{2}^{\prime \prime }+{\it\alpha}^{4}v_{2}].\end{eqnarray}$$
The no-slip and no-penetration conditions at the walls can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle v_{1}={v_{1}}^{\prime }=0\quad \text{at }y=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v_{2}={v_{2}}^{\prime }=0\quad \text{at }y=1. & \displaystyle\end{eqnarray}$$
The kinematic boundary condition gives
$$\begin{eqnarray}h=\frac{v_{1}}{\text{i}{\it\alpha}(U_{1}-c)}=\frac{v_{2}}{\text{i}{\it\alpha}(U_{2}-c)}\quad \text{at }y=h.\end{eqnarray}$$
The continuity of the velocity components across the interface is expressed as
$$\begin{eqnarray}\displaystyle & \displaystyle {v_{1}}^{\prime }-\text{i}{\it\alpha}h{U_{1}}^{\prime }={v_{2}}^{\prime }-\text{i}{\it\alpha}h{U_{2}}^{\prime }\quad \text{at }y=h, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v_{1}=v_{2}\quad \text{at }y=h. & \displaystyle\end{eqnarray}$$
The normal stress jump and continuity of the tangential stress balance in the streamwise and spanwise directions are respectively given by
$$\begin{eqnarray}\displaystyle & & \displaystyle \text{i}{\it\alpha}\mathit{Re}[\{{v_{2}}^{\prime }(c-U_{2})+{U_{2}}^{\prime }v_{2}\}-\{{v_{1}}^{\prime }(c-U_{1})+{U_{1}}^{\prime }v_{1}\}]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\text{e}^{R_{s}}[{v_{2}}^{\prime \prime \prime }-3{\it\alpha}^{2}{v_{2}}^{\prime }]-[v_{1}^{\prime \prime \prime }-3{\it\alpha}^{2}v_{1}^{\prime }]={\it\alpha}^{4}{\it\Gamma}\frac{({v_{2}}^{\prime }-{v_{1}}^{\prime })}{\text{i}{\it\alpha}(U_{2}^{\prime }-U_{1}^{\prime })}\quad \text{at }y=h,\end{eqnarray}$$
$$\begin{eqnarray}\text{e}^{R_{s}}[{v_{2}}^{\prime \prime }+{\it\alpha}^{2}v_{2}]-\frac{(\text{e}^{R_{s}}{U_{2}}^{\prime \prime }-{U_{1}}^{\prime \prime })}{(U_{2}-c)}v_{2}={v_{1}}^{\prime \prime }+{\it\alpha}^{2}v_{1}\quad \text{at }y=h.\end{eqnarray}$$
Here
${\it\Gamma}\equiv {\it\gamma}H/{\it\mu}_{1}Q$
is an inverse capillary number, in which
${\it\gamma}$
denotes the interfacial tension. The complete derivation and linearisation of the stability equations can be found in Sahu & Matar (Reference Sahu and Matar2010). In this work, we set
${\it\Gamma}$
to zero, because we wish to compare the miscible and immiscible cases without the additional factor of surface tension in the latter.
4 Budget of disturbance kinetic energy
An energy budget analysis can highlight the physical differences between the two flows in their stability behaviour. A budget of disturbance kinetic energy, neglecting the surface-tension and gravity, is given by
$$\begin{eqnarray}2{\it\omega}_{i}\frac{1}{{\it\lambda}}\int _{a}^{b}\int _{0}^{{\it\lambda}}E\,\text{d}x\,\text{d}y=\frac{1}{{\it\lambda}}\int _{a}^{b}\int _{0}^{{\it\lambda}}P\,\text{d}x\,\text{d}y+\frac{1}{{\it\lambda}\mathit{Re}}\int _{a}^{b}\int _{0}^{{\it\lambda}}D\,\text{d}x\,\text{d}y+I,\end{eqnarray}$$
where
${\it\lambda}\equiv 2{\rm\pi}/{\it\alpha}$
. For fluid 1 and fluid 2 (
$a=h,b=1$
) and (
$a=0,b=h$
), respectively. The kinetic energy, the rate of its production and the rate of dissipation are given respectively by
$$\begin{eqnarray}E={\displaystyle \frac{1}{2}}(u^{2}+v^{2}),\quad P=-uv\frac{\text{d}U}{\text{d}y},\end{eqnarray}$$
and
$$\begin{eqnarray}D=2{\it\mu}\left[\left(\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial v}{\partial y}\right)^{2}+\frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^{2}\right].\end{eqnarray}$$
The viscous work done by the mean flow on the interface is given by
$$\begin{eqnarray}I=\frac{1}{{\it\lambda}\mathit{Re}}\int _{0}^{{\it\lambda}}\left[u_{1}{\it\tau}_{1}^{xy}-u_{2}{\it\tau}_{2}^{xy}\right]\text{d}x\quad (\text{at }y=h),\end{eqnarray}$$
wherein
$$\begin{eqnarray}{\it\tau}^{xy}={\it\mu}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right).\end{eqnarray}$$
Continuity of shear stresses implies that for the mean flow there is a jump in the slope of
$U$
, i.e.,
$$\begin{eqnarray}U_{1}^{\prime }\big|_{y=h}=\text{exp}(R_{s})U_{2}^{\prime }\big|_{y=h},\end{eqnarray}$$
and for the disturbance
$$\begin{eqnarray}{\it\tau}_{1}^{xy}\big|_{y=h}={\it\tau}_{2}\big|_{y=h}^{yx}={\it\tau}^{xy}.\end{eqnarray}$$
Thus, equation (4.4) can be written as
$$\begin{eqnarray}I=\frac{1}{{\it\lambda}\mathit{Re}}\int _{0}^{{\it\lambda}}{\it\tau}^{xy}\left[u_{1}-u_{2}\right]\text{d}x\quad (\text{at }y=h).\end{eqnarray}$$
When the interface is being deformed, streamwise disturbance velocities of unequal size are forced at the interface, i.e.,
$u_{1}\neq u_{2}$
, due to which energy transfer occurs from the mean flow to the disturbance. This quantity will remain positive even if the fluid layers are interchanged (Boomkamp & Miesen Reference Boomkamp and Miesen1996).
For miscible flow, we would have the same expressions for
$E$
,
$P$
and
$D$
, but integrated across the entire channel, and of course
$I=0$
.
Next we evaluate how miscible and immiscible two-fluid flows differ in their stability behaviour. We then perform direct numerical simulations and show that the nonlinear behaviour is consistent with the predictions of linear instability. The simulations also help us to estimate how three-dimensional the flow is.
5 Results
5.1 Linear stability analysis
A log viscosity ratio of
$R_{s}>1$
gives rise to a velocity profile with a point of inflexion for
$h<0.5$
, as seen in figure 2. Such a profile is likely to be more unstable than one without a point of inflexion, and therefore be the more interesting case, so we restrict ourselves to positive values of
$R_{s}$
. A typical set of disturbance growth rates is presented in figure 3. The growth rates of the most unstable eigenmode are plotted as functions of the wavenumber, for different values of Reynolds number. The instability behaviour for both miscible and immiscible two-fluid flows is shown in the same figure. As is usual in shear flows, the instability gets more severe as the Reynolds number increases. More remarkable is the fact that, for higher Reynolds numbers
$(\mathit{Re}\geqslant 200)$
, the miscible flow is more unstable than the flow containing the immiscible interface.
Figure 3. Growth rates,
${\it\omega}_{i}$
, of the most unstable disturbances as functions of the wavenumber,
${\it\alpha}$
, for different values of the Reynolds number
$\mathit{Re}$
for
$\mathit{Sc}=10$
,
$q=0.02$
,
$R_{s}=1$
and
$h=0.15$
. The solid lines represent miscible flow, and the red dotted lines represent the result for immiscible two-fluid flow, with the interface placed at the same value of
$h$
, and with the same
$R_{s}$
. The symbols for a given Reynolds number are the same for miscible and immiscible flows.
Figure 4. The dispersion curves (
${\it\omega}_{i}$
versus
${\it\alpha}$
) for different values of (a)
$\mathit{Sc}$
for
$q=0.02$
, and (b)
$q$
for
$\mathit{Sc}=10$
. The rest of the parameter values are
$\mathit{Re}=500$
,
$R_{s}=1$
and
$h=0.15$
. The dotted lines represent the results for the immiscible case.
Figure 5. Dispersion curves (
${\it\omega}_{i}$
versus
${\it\alpha}$
) for different values of (a)
$h$
for
$R_{s}=1$
, (b)
$R_{s}$
for
$h=0.3$
. The rest of the parameter values are
$\mathit{Re}=500$
,
$\mathit{Sc}=100$
and
$q=0.01$
. The red dotted lines show the corresponding growth rates in the case of immiscible flow.
It is now accepted knowledge that shear flows of two or more fluids most often become more unstable at high Schmidt numbers (when the diffusivity of one fluid in another is very low). The expectation therefore would be that if we increase the Schmidt number of the miscible flow, flow would become increasingly unstable. We see in figure 4(a) that the behaviour is not monotonic with increased Schmidt number. While the flow becomes more unstable as we increase the Schmidt number up to a value of
$100$
, a further increase in
$\mathit{Sc}$
decreases the growth rate of the most unstable mode. For very high
$\mathit{Sc}$
, i.e. for
$\mathit{Sc}>10^{5}$
, the behaviour of the most unstable mode is the same as that of the immiscible flow. Thus the immiscible case is less unstable than two-fluid flow of intermediate miscibility. We will show later in this section that the overlap of the mixed layer with the critical layer is the underlying mechanism in the present system, which is characteristic of high Reynolds number flow. A non-monotonic response to changes in the Schmidt number was also obtained by Ern et al. (Reference Ern, Charru and Luchini2003) in Couette flow at low Reynolds number, but their mechanism was not the same as this one, as will be discussed below.
Another parameter which is known to affect flow stability significantly is the thickness
$q$
of the mixed layer. In figure 4(b), we see that as the mixed layer is made thinner, the growth rate of the dominant instability increases. This is as expected, and is caused by the fact that as
$q$
decreases, the viscosity gradient becomes sharper, making the stability operator more singular. In this figure,
$\mathit{Sc}=10$
is used as a typical example. It can be seen that the growth rate remains sensitive to
$q$
at all values of
$q$
that we have considered. Here too it can be observed that the dispersion curve for the immiscible flow (shown by the dashed line) is well below the dispersion curves of the miscible system for
$q\leqslant 0.05$
. This figure is for a Reynolds number of
$500$
, but we have repeated all our calculations at a Reynolds number of
$1000$
as well (not shown), and the behaviour is qualitatively the same. Again, when the layer is thin enough, the flow of intermediate miscibility is significantly more unstable than the immiscible case.
In figure 5(a), we investigate the effect of
$h$
, the height of the mixed layer from the bottom wall. When the interfacial layer is close to the bottom or top walls, the flow is stabler than when the mixed layer is near the middle, and a value of
$h\sim 0.3$
is the least stable. The response to the location of the interfacial layer is thus non-monotonic. For
$h<0.4$
we see that the miscible flow is more unstable than the immiscible. It can be seen in figure 5(b) that immiscible flow is not very sensitive to viscosity ratio, but disturbances in miscible flow grow much faster at higher viscosity ratios. For all the viscosity ratios considered, it is seen that the miscible flow is more unstable than the corresponding immiscible flow. Taking into consideration all the linear stability results, we see that our finding that the miscible flow (at intermediate levels of miscibility) is more unstable than the immiscible flow is a general result for high Reynolds number channel flow of two fluids.
In order to investigate the instability mechanism, neutral stability curves for different values of
$\mathit{Sc}$
and
$h$
are plotted in figures 6(a) and 6(b), respectively. The standard Tollmien–Schlichting mode contributes a region of instability, seen on the extreme right of the plots, i.e., at high Reynolds number. In addition, a distinct region of instability is observed, which grows in size with increasing Schmidt number (figure 6
a). The phase speed in this regime is close to the mean velocity in the mixed-fluid layer (we shall return to this point in figure 7). Note that the neighbourhood of thickness
$O(R^{-1/3})$
, where the phase speed of the dominant disturbance is close to the mean velocity, is the critical layer where most of the disturbance kinetic energy is produced. It was shown in Govindarajan (Reference Govindarajan2004) for a three-layer channel flow that the above condition, of an overlap between the critical layer with the mixed layer, contributes to a singular perturbation term in the stability operator. The resulting new mode of instability was termed the ‘overlap’ mode. The energy production is interfered with in a major way by this overlap. Since this instability is inherently inertial it is distinct from the modes obtained by Talon & Meiburg (Reference Talon and Meiburg2011) for Stokes flow. Besides the fact that Talon & Meiburg found a similarity between their instability and that of Ern et al. (Reference Ern, Charru and Luchini2003), we may check directly whether there is an overlap mechanism operational in the latter. It can be checked that the critical layer obtained by Ern et al. (Reference Ern, Charru and Luchini2003) is well below their mixed layer. Thus the instabilities observed by Talon & Meiburg (Reference Talon and Meiburg2011) and Ern et al. (Reference Ern, Charru and Luchini2003) are not overlap modes. In fact in Stokes flow, the dissipation of the overlap mode would be infinite and energy production could never exceed dissipation in order to make the flow unstable.
Figure 6. Neutral stability curves for different values of (a)
$\mathit{Sc}$
for
$h=0.8$
, and (b)
$h$
for
$\mathit{Sc}=1000$
. The rest of the parameter values are
$R_{s}=0.6$
and
$q=0.05$
. For the boundaries which appear as closed curves, the region contained within is unstable, and the outside is stable. For the open curves, the region to the right is unstable, and that to the left is stable. The curves at high Reynolds number (of approximately 10 000, where
${\it\alpha}\sim 2$
) correspond to the Tollmien–Schlichting mode of instability, whereas the other curves correspond to the overlap mode, as defined in Govindarajan (Reference Govindarajan2004).
It is seen in figure 6(b) that the overlap mode displays a distinct region of instability when
$h\geqslant 0.7$
, whereas for lower values of
$h$
a much larger region is unstable, and it is difficult to distinguish the ‘overlap’ mode anymore, except that it can be recognised by an apparent kink in the neutral boundary. This behaviour is observed over a range of parameters, and an example at low
$\mathit{Sc}$
is shown in figure 7(a) for
$R_{s}=1$
,
$h=0.2$
and
$q=0.05$
. To verify whether some part of the neutral boundary corresponds to an overlap mode of instability, we examine figure 7(b), which represents behaviour along the lower limb of the neutral stability boundary of figure 7(a). The distance between the centres of the mixed layer and the critical layer is the quantity
$h-y_{cr}$
. This quantity is plotted versus Reynolds number along the lower limb of the neutral stability boundary in this figure. The width of the critical layer may be estimated as
${\sim}(U_{cr}^{\prime }\mathit{Re}{\it\alpha})^{-1/3}$
, and this is denoted by the region within the red lines. It is now evident that different modes of instability are in operation on either side of the kink seen in the neutral stability boundary in figure 7(a). The neutral mode at low Reynolds numbers has a small distance between
$y_{cr}$
and
$h$
. In fact this distance is seen to be smaller than the critical layer thickness over a range of Reynolds numbers, indicating that overlap effects must be in operation at the lowest order. At high Reynolds number however (beyond the kink, coming downwards along the neutral boundary), a sudden jump is seen in
$h-y_{cr}$
, and this difference is greater than the critical layer thickness, indicating that different physics is operational there.
Figure 7. (a) Neutral stability curves, (b)
$y_{cr}-h$
versus
$\mathit{Re}$
along the lower limb of the neutral stability curves (shown in (a)). Without symbols:
$\mathit{Sc}=0.01$
, with symbols:
$\mathit{Sc}=1$
. The rest of the parameter values are
$R_{s}=1$
,
$h=0.2$
and
$q=0.05$
. The regions contained within the red lines correspond to an order of magnitude estimate of the critical layer thickness
$\pm (U_{cr}^{\prime }\mathit{Re}{\it\alpha})^{-1/3}/2$
. Here
$y_{cr}$
is the location of the critical layer, i.e. the
$y$
at which
$U=c_{r}$
.
The production of disturbance kinetic energy is examined next for the intermediate Schmidt number case of figure 3. The variations of the disturbance kinetic energy rate
$E$
, the production rate
$P$
and the dissipation rate
$D$
across the channel are displayed in figures 8(a), 8(b) and 8(c), respectively. At
$R_{s}=1$
, the maximum growth rate in both miscible and immiscible flows occurs at a wavenumber
${\it\alpha}\sim 4.5$
, so the most dangerous mode at
${\it\alpha}=4.5$
is chosen to do this energy budget analysis. The kinetic energy production is seen in figure 8(c) to peak close to the location of the mixed layer, indicating that overlap effects are in operation. The difference between the production and the dissipation rates, which gives, in the miscible flow case, the change disturbance kinetic energy per unit time, is shown in figure 8(d). The production and the dissipation in the portions of the channel away from the interfacial or mixed layer are very similar in the two flows. The major difference is apparent in the vicinity of the interfacial or the viscosity stratified layer. It is clear that the net production
$P-D$
in the miscible case exceeds that in the immiscible case in the mixed region. The immiscible flow has an additional contribution
$I$
to disturbance growth. We denote by
${\mathcal{E}}$
,
${\mathcal{P}}$
and
$\mathscr{D}$
the integrals of
$E$
,
$P$
and
$D$
across the channel from wall to wall. It is seen from (4.1) that the growth rate is given by
$2{\it\omega}_{i}=({\mathcal{P}}-\mathscr{D}+I)/{\mathcal{E}}$
.
${\mathcal{P}}$
,
${\mathcal{D}}$
and
$I$
normalised with
${\mathcal{E}}$
for the immiscible case are 0.7479,
$-0.6054$
and 1.03423, respectively. This gives the growth rate
${\it\omega}_{i,max}=0.5878$
, where the subscript
$max$
stands for the maximum growth rate of the fastest-growing mode. For the miscible case: the values of
${\mathcal{P}}$
and
${\mathcal{D}}$
normalised with
${\mathcal{E}}$
are 5.6366 and
$-3.5407$
, respectively, giving
${\it\omega}_{i,max}=1.045$
. It is thus seen that the net production minus dissipation on disturbance kinetic energy in the miscible case is larger than the contribution of all terms in the immiscible case.
Figure 8. Variation of (a) kinetic energy
$E$
, (b) dissipation rate. The negative (
$-D$
) is shown to help viewing. (c) Production rate
$P$
and (d)
$P-D$
of the most dangerous disturbance mode (
${\it\alpha}=4$
) in the wall-normal direction for
$\mathit{Re}=500$
. The rest of the parameters are the same as those used to generate figure 3. The values of
${\mathcal{P}}$
,
${\mathcal{D}}$
and
$I$
normalised with
${\mathcal{E}}$
for the immiscible case are 0.8458,
$-1.0745$
and 1.1326, respectively;
${\it\omega}_{i,max}=0.2159$
. For the miscible case, the values of
${\mathcal{P}}$
and
$\mathscr{D}$
normalised with
${\mathcal{E}}$
are 4.5379 and
$-1.9227$
, respectively;
${\it\omega}_{i,max}=0.3269$
.
5.2 Three-dimensional numerical simulations
5.2.1 Numerical method
For miscible systems, equations (2.3)–(2.5) are solved by a finite-volume approach (Ding, Spelt & Shu Reference Ding, Spelt and Shu2007) using a staggered grid discretisation, i.e. the scalar variables (the pressure and concentration of the scalar) and the velocity components are defined at the centre and at the cell faces, respectively. The discretised convection-diffusion equation of
$s_{0}$
is given by
$$\begin{eqnarray}\frac{\frac{3}{2}{s_{0}}^{n+1}-2{s_{0}}^{n}+\frac{1}{2}{s_{0}}^{n-1}}{{\rm\Delta}t}=\frac{1}{\mathit{Re}\mathit{Sc}}{\rm\nabla}^{2}{s_{0}}^{n+1}-2\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}^{n}{s_{0}}^{n})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}^{n-1}{s_{0}}^{n-1}),\end{eqnarray}$$
where
${\rm\Delta}t=t^{n+1}-t^{n}$
and the superscript
$n$
represents the time step. The advective terms, i.e. the nonlinear terms in (2.5) are discretised using a weighted essentially non-oscillatory scheme, and a central difference scheme is used to discretise the diffusive terms on the right-hand side of (2.4)–(2.5). Second-order accuracy in the temporal discretisation is obtained by employing the Adams–Bashforth and the Crank–Nicolson methods for the advective and second-order dissipation terms in (2.4), respectively. This discretised form of (2.4) is given by
$$\begin{eqnarray}\displaystyle \frac{\boldsymbol{u}^{\ast }-\boldsymbol{u}^{n}}{{\rm\Delta}t}=\frac{1}{p^{n+1/2}}\left\{-\left[\frac{3}{2}{\mathcal{H}}(\boldsymbol{u}^{n})-\frac{1}{2}{\mathcal{H}}(\boldsymbol{u}^{n-1})\right]+\frac{1}{2\mathit{Re}}\left[{\mathcal{L}}(\boldsymbol{u}^{\ast },{\it\mu}^{n+1})+{\mathcal{L}}(\boldsymbol{u}^{n},{\it\mu}^{n})\right]\right\},\hspace{-15.00002pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}^{\ast }$
is the intermediate velocity, and
${\mathcal{H}}$
and
${\mathcal{L}}$
denote the discrete convection and diffusion operators, respectively. The intermediate velocity
$\boldsymbol{u}^{\ast }$
is then corrected to
$(n+1)\text{th}$
time level.
$$\begin{eqnarray}\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{\ast }}{{\rm\Delta}t}=\boldsymbol{{\rm\nabla}}p^{n+1/2}.\end{eqnarray}$$
The pressure distribution is obtained from the continuity equation at time step
$n+1$
using
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\boldsymbol{{\rm\nabla}}p^{n+1/2}\right)=\frac{\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }}{{\rm\Delta}t}.\end{eqnarray}$$
For the corresponding immiscible system, a diffuse interface method based on Ding et al. (Reference Ding, Spelt and Shu2007) is used. In this case, instead of the diffusion equation (i.e. (2.5)), the Cahn–Hilliard equation, given by
$$\begin{eqnarray}\frac{\partial s_{0}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}s_{0}=\frac{1}{\mathit{Re}\mathit{Sc}_{i}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(M\boldsymbol{{\rm\nabla}}{\it\phi}),\end{eqnarray}$$
where
$\mathit{Sc}_{i}\equiv {\it\mu}_{1}/({\it\rho}M_{c}{\it\phi}_{c})$
, is solved. Here
$M_{c}$
and
${\it\phi}_{c}$
are the characteristic values of mobility and chemical potential,
${\it\phi}$
(
$\equiv {\it\epsilon}^{-1}{\it\sigma}{\it\alpha}{\it\Psi}^{\prime }(s_{0})-{\it\epsilon}{\it\sigma}{\it\alpha}{\rm\nabla}^{2}s_{0}$
), respectively, wherein
${\it\epsilon}$
is the measure of interface thickness,
${\it\Psi}(s_{0})={s_{0}}^{2}(1-s_{0})^{2}/4$
is the bulk energy density and
${\it\alpha}$
is a constant.
The above discretised equations are solved by employing the no-slip and the no-penetration boundary conditions for the velocity components and no-flux condition for concentration
$(s_{0})$
at the walls, and periodic boundary conditions in the axial and lateral directions. The pressure gradient is kept the same as that of the stability analysis conducted in the previous section. A domain with 321, 81 and 161 cells in the axial
$x$
, wall-normal
$y$
and spanwise
$z$
, directions, respectively, are used in the simulations. Grid refinement tests have been conducted to ensure that we obtain grid-converged results. The numerical procedures described above for miscible and immiscible systems are similar to the ones used by Ding et al. (Reference Ding, Spelt and Shu2007). The reader is referred to this paper for detailed descriptions and validation of the solvers.
In figure 9, we compare the maximum value of the wall-normal velocity component,
$v_{max}$
, obtained from our direct numerical simulations with that obtained from linear stability analysis for
$\mathit{Re}=500$
,
$\mathit{Sc}=100$
,
$R_{s}=2$
,
$h=0.3$
and
$q=0.01$
. The miscible and immiscible simulations are conducted inside a channel with a length that is twice the wavelength of the most dangerous mode obtained from linear stability analysis (
${\it\alpha}=4$
for this set of parameter values). If we were to prescribe the dominant perturbation mode in the initial conditions, we would expect this perturbation to grow in consonance with linear theory at early times. We choose to study the harder case, i.e. one in which we prescribe no initial perturbation. Thus all modes of perturbation are equally initialised, so the dominant one will need time to become visible and to grow in accordance with linear theory. In the miscible case, we see that at some time (after
$t=1$
), the growth rate of the disturbance in the numerical simulations is close to that of the dominant mode predicted by linear theory (seen by the fact that the two lines are parallel). In the immiscible case, the dominant mode is not distinguishable by its linear growth rate at any time, indicating that nonlinear effects dominate the entire process. In the nonlinear regime too, it can be observed that the miscible flow with finite Schmidt number is more unstable than the immiscible flow. The subsequent three-dimensional plots (figures 10 and 11) convey the same information but in pictorial form of an iso-concentration contour.
Figure 9. Comparison of the temporal evolution of the vertical perturbation obtained from linear stability analysis (solid lines) and direct numerical simulations (dashed lines). Miscible flow: symbols; immiscible flow: no symbols. The rest of the parameter values are
$\mathit{Re}=500$
,
$\mathit{Sc}=100$
,
$R_{s}=2$
,
$h=0.3$
and
$q=0.01$
.
Figure 10. Temporal evolution of (a,c,d,g,h)
$s=0.5$
contour for miscible flow with
$q=0.05$
,
$\mathit{Sc}=100$
; (b,e,f,i,j) interface for the immiscible case. (a,b)
$t=4$
, (c–f) 6 and (g–j) 10. The rest of the parameter values are
$h=0.3$
,
$\mathit{Re}=500$
and
$R_{s}=2$
. The flow is along the positive
$x$
direction. The side panels at
$t=6$
are the cross-sectional views in the
$x$
–
$z$
plane at
$y=1$
, and those at
$t=10$
are the cross-sectional views in the
$y$
–
$z$
plane at
$x=3$
.
Figure 11. Temporal evolution of (a,c,e,f,i)
$s=0.5$
contour for miscible flow with
$q=0.01$
,
$\mathit{Sc}=100$
; (b,d,g,h,j) interface for the immiscible case. (a,b)
$t=3$
, (c,d) 5, (e–h) 7 and (i,j) 9. The rest of the parameter values are
$h=0.4$
,
$\mathit{Re}=500$
and
$R_{s}=2$
. The flow is along the positive
$x$
direction. The insets at
$t=7$
represent cross-sectional views in the
$x$
–
$z$
plane.
A few cases were computed, of which we present results for a Reynolds number of
$\mathit{Re}=500$
and
$R_{s}=2$
as being representative. The spatio-temporal evolution of the interface separating the fluids in the immiscible flow (
$q=0$
,
$\mathit{Sc}=\infty$
) and the
$s=0.5$
contour for the miscible flow (
$\mathit{Sc}=100$
) for
$q=0.01$
,
$h=0.4$
and
$q=0.05$
,
$h=0.3$
are presented in figures 10 and 11, respectively. A computational domain of size
$4\times 2\times 1$
is used for these simulations, wherein velocity components are set to zero initially. Given the constraints of numerical accuracy, we did not perform simulations with thinner mixed layers. It can be seen in figure 10(a,c,d,g,h) that the miscible flow becomes unstable, and rolled-up structures are obtained, which can be observed at
$t=6$
in the contour of
$s=0.5$
. At later times the nonlinear instability develops an irregular lateral structure too (see the side panel at
$t=10$
in figure 10
g,h). However, it can be seen in figure 10(b,e,f,i,j) that the interfacial flow becomes unstable at a time much later than the corresponding miscible flow. Secondly, the immiscible flow remains two-dimensional in the regime where the miscible flow has become three-dimensional. For
$h=0.4$
(figure 11) the contrast is even more pronounced; in this case the interfacial flow is stable till
$t=9$
, whereas the miscible flow becomes unstable at
$t\approx 5$
. At early times this behaviour is consistent with that obtained in the linear stability analysis, and at later times it is seen that the tendency for the immiscible flow to remain less unstable than the miscible persists into the nonlinear regime as well.
6 Discussion and summary
In the present study, we have investigated the difference between interfacial flow instabilities (without surface tension) with a viscosity jump across the interface, and instabilities associated with two-layer miscible channel flow of two fluids with different viscosities. We show that in this flow too, an overlap mode of instability, seen before in other miscible flow configurations (Govindarajan & Sahu Reference Govindarajan and Sahu2014) is dominant, where the mixed layer overlaps significantly with the critical layer of the dominant disturbance. This means that the lower-order terms in the critical layer balance get disturbed by the viscosity variation in the mixed layer and contribute to significant changes in the stability behaviour. The overlap mode of instability is produced by an inertial effect, and is driven by different physics from instabilities seen before in this flow under zero Reynolds number conditions. In fact the overlap mode of instability will necessarily vanish at zero Reynolds number. At Reynolds numbers of order
$1$
or lower, the critical layer is as wide as the flow, and is thus unable to produce singular effects. Interestingly the overlap mode of instability becomes operational at relatively low Schmidt numbers and low viscosity ratios, unlike the mechanisms in operation at very low Reynolds number, which are triggered by high viscosity ratio and poor diffusivity.
We have studied how decreasing the diffusivity of the fluids takes the results closer to the immiscible case. Above a Schmidt number of
${\sim}10^{5}$
, the behaviour of the miscible layer is very close to that of the immiscible. At moderate Schmidt number, of up to
$100$
, we find that the channel flow of two miscible fluids can be more unstable than the case where the two fluids are immiscible. The increase in the disturbance growth rate due to finite diffusivity is not too large, but makes a point of principle which needs further investigation. Whether this effect at high Reynolds number has any connection with the destabilisation due to diffusivity seen at very low Reynolds numbers by Ern et al. (Reference Ern, Charru and Luchini2003) is not established yet, but as discussed above we have strong reasons to believe that the two are completely different.
As expected, increasing the Reynolds number destabilises the flow. We find that increasing the viscosity contrast between the two fluids does not have any significant effect on instability characteristics in the immiscible case, but significantly increases the growth rate for the miscible two-layer flow. Reducing the thickness of the mixed layer increases the growth rate, as expected. Thus for thin mixed layers at intermediate diffusivity, the increase in instability as compared to the immiscible case is larger. Varying the location of the mixed layer or the interface has an effect on the stability. In some range of this parameter, distinct regions of overlap instability are obtained, whereas in others, the instability regions due to different mechanisms merge with each other.
The miscible two-layer channel flow has not been studied earlier in the nonlinear regime to our knowledge, and we therefore conduct direct numerical simulations for both the miscible and immiscible cases. At a Schmidt number of 100, linear stability analysis predicts a faster growth rate for the miscible than for the immiscible. This is borne out by the simulations. Also, we see that nonlinear effects on the immiscible flow are visible at earlier times than in the miscible, and the rate at which we see the interface roll-up can be made much slower by immiscibility, or even suppressed.
Acknowledgement
The authors would like to sincerely thank the anonymous reviewers for their valuable comments and suggestions.
Appendix A. Validity of the parallel flow assumption
Consider a situation in which a splitter plate is located at
$x<x_{0}$
, at a constant
$y$
, and parallel streams of two miscible fluids flow on both sides of this plate. The streams come into contact with each other at
$x=x_{0}$
. The two fluids begin to mix with each other for
$x>x_{0}$
, thus producing a stratified layer. The thickness,
$q$
, of this layer grows as the fluids move in the downstream direction and therefore
$q$
is a function of
$x$
. We note that the flow diffuses as it moves downstream but does not diffuse in time at one
$x$
location, i.e., the base flow is steady in time.
We know that at any location, the concentration
$s_{0}$
satisfies the following equation,
$$\begin{eqnarray}\frac{\partial s_{0}}{\partial t}+U\frac{\partial s_{0}}{\partial x}+V\frac{\partial s_{0}}{\partial y}=\frac{1}{\mathit{Re}\mathit{Sc}}\left[\frac{\partial ^{2}s_{0}}{\partial x^{2}}+\frac{\partial ^{2}s_{0}}{\partial y^{2}}\right].\end{eqnarray}$$
For slow diffusion (i.e for high
$Pe\equiv \mathit{Re}\mathit{Sc}$
), we can make the assumption on locally parallel flow (variation of
$s_{0}$
in the
$y$
direction is much larger than that in the
$x$
direction); thus
$V\ll U$
and
$\partial ^{2}/\partial x^{2}\ll \partial ^{2}/\partial y^{2}$
. This is equivalent to saying that the variations of the gradients of the flow variables and the thickness
$q$
of the mixed region have much larger length scale than the disturbance wavelength. In such a scenario, the concentration is a function of
$y$
and
$t$
only and not of
$x$
. Thus the above equation reduces to
$$\begin{eqnarray}\frac{\partial s_{0}}{\partial t}=\frac{1}{Pe}\frac{\partial ^{2}s_{0}}{\partial y^{2}}.\end{eqnarray}$$
Using the same approximation, we know that
$U\sim O(1)$
,
$y\sim \sqrt{{\it\nu}}$
, where
${\it\nu}$
is the kinematic viscosity, as the viscosity is directly proportional to the concentration in the mixed layer. Therefore,
$qs_{0}\sim O(y^{2})$
. This implies that
$\partial s_{0}/\partial x\simeq (1/q)\,O(1/Pe)$
. Thus for
$Pe\gg 1$
,
$\partial s_{0}/\partial x$
is very small, i.e., the downstream variation of
$s_{0}$
is very small which in turn implies that the change in the thickness of the mixed layer
$(q)$
along the
$x$
-direction is very small.
Alternatively, if we assume a similarity solution
$s_{0}(y/q(x))\simeq s({\it\xi})$
(where
${\it\xi}=(y/q(x))$
), from (A 1), we get
$$\begin{eqnarray}U\frac{\text{d}s_{0}}{\text{d}{\it\xi}}\left(-\frac{{\it\xi}}{q}\,\frac{\text{d}q}{\text{d}x}\right)\simeq \frac{1}{Pe}\left(\frac{\text{d}^{2}s_{0}}{\text{d}{\it\xi}^{2}}\,\frac{1}{q^{2}}\right).\end{eqnarray}$$
As a consequence,
$$\begin{eqnarray}\frac{1}{q}\,\frac{\text{d}q}{\text{d}x}\sim \frac{1}{q^{2}Pe}\Rightarrow \frac{\text{d}q}{\text{d}x}\sim \frac{1}{q}\,O(Pe)^{-1}.\end{eqnarray}$$
Thus, the downstream growth of the mixed layer is inversely proportional to the Péclet number as
$U$
and
${\it\xi}$
are of
$O(1)$
, and
$O(\text{d}s_{0}/\text{d}{\it\xi})\simeq O(\text{d}^{2}s_{0}/\text{d}{\it\xi}^{2})$
. For most of the simulations considered in the present study,
$Pe\geqslant 1000$
,
$\mathit{Re}\geqslant 100$
and
$q\geqslant 0.05$
; these parameter values are well above the limit for which parallel flow assumption is valid. This confirms that for the Reynolds and the Schmidt numbers considered in the present study, the assumption of uniform thickness of the viscosity stratified layer is justified.
Figure 12. (a)
$S_{0}$
profile, (b) the dispersion curve (
${\it\omega}_{i}$
versus
${\it\alpha}$
) obtained using the present base state and an error function profile. The other parameters are chosen as
$\mathit{Re}=500$
,
$\mathit{Sc}=10$
,
$R_{s}=1$
,
$h=0.3$
and
$q=0.1$
.
Now, the solution of (A 2) is an error function, and the fifth-order polynomial is a good representation of this, as seen in figure 12(a) in this response. The growth rates obtained using an error function type profile and the fifth-order polynomial used in the present study are compared in figure 12(b). It can be seen that they too agree very well.
We also note that such basic flows are commonly used in stability studies. Several authors have used the same logic to give a basic concentration profile in the form of a hyperbolic tangent (Ern et al. Reference Ern, Charru and Luchini2003) or an error function (e.g. Selvam et al. Reference Selvam, Talon, Lesshafft and Meiburg2009; Talon & Meiburg Reference Talon and Meiburg2011), given by
$$\begin{eqnarray}s_{0}=0.5-0.5\text{erf}\left[\frac{y-h-0.5q}{q}\right].\end{eqnarray}$$
Some have also used a fifth-order polynomial (see e.g. Malik & Hooper Reference Malik and Hooper2005) which is smooth enough to approximate either profile.
In summary, the disturbance wavelength is much shorter than the downstream length scale over which the mixed-layer thickness registers any growth, so it is justifiable to use locally a constant-thickness approximation.
