1 Introduction
The stability of plane Couette–Poiseuille two-fluid flows was first investigated by Yih (Reference Yih1967), who showed that the flow is unstable to long waves at non-zero Reynolds numbers. This instability is due to a viscosity jump across the interface, and depends on the viscosity and thickness ratios of the two fluid layers (the density and surface tension ratios were held fixed in Yih’s study). Further work by Hooper & Boyd (Reference Hooper and Boyd1983), Hooper (Reference Hooper1985) and Renardy (Reference Renardy1985, Reference Renardy1987) revealed that the flow is always linearly unstable if the thin layer is also more viscous; this effect is known as the ‘thin-layer effect’ and is often found in multilayer flows (Craster & Matar Reference Craster and Matar2009).
The aim of the present study is the analysis of the nonlinear stages of these instabilities and a direct comparison of a new theory with both direct numerical simulations (DNS) and the experiments of Barthelet, Charru & Fabre (Reference Barthelet, Charru and Fabre1995). Previous attempts have been made using Kuramoto–Sivashinsky (KS) equations – see for example Hooper & Grimshaw (Reference Hooper and Grimshaw1985, Reference Hooper and Grimshaw1988) and subsequently Charru & Fabre (Reference Charru and Fabre1994) and Barthelet & Charru (Reference Barthelet, Charru, Renardy, Coward, Papageorgiou and Sun1995), who take smaller surface tension and retain higher-order terms corresponding to linear and nonlinear dispersion as well as nonlinear dissipation. Barthelet & Charru (Reference Barthelet, Charru, Renardy, Coward, Papageorgiou and Sun1995) present a detailed and critical evaluation of KS-type equations in the light of the experiments by the same authors. The experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995) are set up in an annular Couette device which facilitates the generation of periodic waves (described in more detail in § 2), and predict that as the upper plate speed increases, and for thickness ratios less than approximately 0.3 (the thin layer is always more viscous), there is a supercritical bifurcation to non-symmetric nonlinear travelling waves with wavelength equal to the Couette device length. In addition, bistability is found at larger plate speeds, with unimodal and bimodal (i.e. wavelengths equal to half the device length) solutions coexisting and emerging as long-time features in the experiments depending on the initial conditions. For order-one thickness ratios, subcritical bifurcations are found with travelling waves of shorter wavelength scaling with the device thickness; this branch of solutions also supports hysteresis, with the travelling wave persisting as the plate velocity is gradually decreased (see, for example, figure 10 in Barthelet et al. (Reference Barthelet, Charru and Fabre1995)). The main focus of the experiments was on the supercritical waves, and detailed data are provided for a thickness ratio of 0.25. We will concentrate on this relatively small thickness ratio regime and develop a novel set of equations capable of describing these nonlinear phenomena. It is worth noting that the detailed evaluation of KS-type models by Charru & Fabre (Reference Charru and Fabre1994) identified some shortcomings of the models vis a vis the experiments: the KS equation is found to be inappropriate (no linear dispersion is present and reflectionally symmetric profiles emerge), while a judicious selection of higher-order terms can reproduce the experimental features, at least qualitatively. The basis of the weakly nonlinear theory of Hooper & Grimshaw (Reference Hooper and Grimshaw1985) and Charru & Fabre (Reference Charru and Fabre1994) is that given a wave amplitude
$\unicode[STIX]{x1D6FC}\ll 1$
, the wavelength is of order
$1/\unicode[STIX]{x1D6FC}\gg 1$
. This introduces an arbitrarily large wavelength and as a result inertial effects cannot be accounted for. Our approach is very different – the wavelength is set by the geometry of the problem (e.g. the channel thickness), and analytical progress is possible when one of the layers is thin. For the thin layer the wave is long but for the thicker one it is not – in fact it scales with the thickness of the respective layer. Asymptotic solutions are sought in each layer and matched at the interface to provide a single evolution equation for the scaled amplitude which accounts for the flow in the thick layer at arbitrary Reynolds numbers and in particular for those of the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995). In addition to non-symmetric profiles, our model also predicts the experimentally observed bistability phenomena. Similar analyses were carried out for pressure-driven two-phase core–annular flows by Papageorgiou, Maldarelli & Rumschitzki (Reference Papageorgiou, Maldarelli and Rumschitzki1990) and by Kalogirou & Papageorgiou (Reference Kalogirou and Papageorgiou2016) for a small variation of the present problem but in the presence of insoluble surfactants. Related non-local equations are found in other physical contexts, for example turbulent shear flows in riverbed dynamics (Fowler Reference Fowler2011), the dynamics of electrified falling films (Tseluiko & Papageorgiou Reference Tseluiko and Papageorgiou2010) or the dynamics of falling films below turbulent gas flows (Tseluiko & Kalliadasis Reference Tseluiko and Kalliadasis2011).
The structure of the rest of this paper is as follows. Section 2 provides a brief overview of the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995) and identifies the ones that are used for direct comparisons. Section 3 presents the asymptotic model and a brief derivation of the weakly nonlinear interfacial evolution equation. In § 4, numerical results obtained by the asymptotic model are compared with DNS and the experiments. Bistability results are also presented. A discussion can be found in § 5.
2 The experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995)
Barthelet et al. (Reference Barthelet, Charru and Fabre1995) observed long-wave instability in a channel of rectangular cross-section bent into an annular ring (see their figure 2). The flow was driven by the rotation of the upper plate, and the position of the interface was measured with a probe situated in the middle of the cross-section. The lower experimental fluid selected here was a mixture of distilled water and glycerine in proportions 42 %–58 % – other mixtures were also examined in the experiments. The upper fluid was a mineral oil. The channel dimensions and the physical properties of the fluids are given in table 1. The channel width to diameter ratio was
$1/10$
and so the two-dimensional flow configuration of figure 1 is appropriate. The experiments were performed for various values of the depth ratio
$h=h_{2}/h_{1}$
, but the case that is most relevant in comparisons with our asymptotic theory is
$h=0.25$
, which is the smallest
$h$
they investigated.
Table 1. Physical parameters from the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995). These are used in the computations of the model equation and the DNS.
The experiments show that for
$h=0.25$
there is a critical upper plate speed above which a supercritical bifurcation occurs to give a travelling wave whose wavelength is equal to the mean channel perimeter
$L_{w}=\unicode[STIX]{x03C0}D=1.257$
m, implying a dimensionless period
$2L=L_{w}/d\approx 63$
. In this paper, we focus on data from three experiments, as detailed in table 2, for three different upper plate speeds
$U$
(the definitions of the Reynolds numbers
$\mathit{Re}_{1}$
,
$\mathit{Re}_{2}$
, the Weber number
$We$
and the capillary number
$Ca$
are given in § 3). As mentioned in the introduction, subcritical shorter waves that scale with
$d$
have also been observed, but at values of the depth ratio
$h$
that are of order one. We do not pursue these here because they are outside the range of validity of the model equations derived below; DNS may be capable of describing them, but this is left for future work.
Figure 1. Geometry of the problem: two superposed fluid layers in a channel of depth
$d$
, driven by the upper plate motion with speed
$U$
.
Table 2. Dimensionless parameters for three different upper plate speeds
$U$
used in the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995) – other parameters are as in table 1.
3 Nonlinear evolution equations in the thin-layer limit
Consider the flow of two immiscible, incompressible and viscous fluids of equal densities
$\unicode[STIX]{x1D70C}_{1}=\unicode[STIX]{x1D70C}_{2}=\unicode[STIX]{x1D70C}$
and different viscosities
$\unicode[STIX]{x1D707}_{1}$
and
$\unicode[STIX]{x1D707}_{2}$
in a channel of height
$d$
, driven by the motion of the upper plate with velocity
$U$
(figure 1). Using Cartesian coordinates, the channel walls are at
$y=0$
and
$y=d$
, and the interface separating the fluids is at
$y=S(x,t)$
, where
$t$
is time. Surface tension is present at the interface and has constant value
$\unicode[STIX]{x1D6FE}$
. In the undisturbed state, the interface is at
$y=\ell d$
, where
$0<\ell <1$
is a constant. The domains
$0<y<S(x,t)$
and
$S(x,t)<y<d$
are denoted as regions 1 and 2 respectively. Lengths are non-dimensionalised with
$d$
, velocities with
$U$
, time with
$d/U$
and pressures with
$\unicode[STIX]{x1D70C}U^{2}$
. The Navier–Stokes and continuity equations in each phase
$i=1,2$
become (gravity is neglected here)
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{i}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{i}=-\unicode[STIX]{x1D735}p_{i}+\frac{1}{\mathit{Re}_{i}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{i},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{i}=0,\end{eqnarray}$$
where
$\boldsymbol{u}_{i}=(u_{i},v_{i})^{\text{T}}$
is the velocity field,
$p_{i}$
is the pressure and
$\mathit{Re}_{i}$
is the Reynolds number in each fluid,
$$\begin{eqnarray}\mathit{Re}_{1}=\frac{\unicode[STIX]{x1D70C}Ud}{\unicode[STIX]{x1D707}_{1}},\quad \mathit{Re}_{2}=\frac{\unicode[STIX]{x1D70C}Ud}{\unicode[STIX]{x1D707}_{2}}=\frac{\mathit{Re}_{1}}{m},\end{eqnarray}$$
with
$m=\unicode[STIX]{x1D707}_{2}/\unicode[STIX]{x1D707}_{1}$
the viscosity ratio. The boundary conditions are no-slip at the walls,
$$\begin{eqnarray}u_{1}=v_{1}=0\quad \text{at}~y=0,\quad \text{and}\quad u_{2}=1,\quad v_{2}=0\quad \text{at}~y=1,\end{eqnarray}$$
along with continuity of velocities at the interface,
$$\begin{eqnarray}u_{1}=u_{2},\quad v_{1}=v_{2}\quad \text{at}~y=S(x,t).\end{eqnarray}$$
In addition, at the interface we need to satisfy a kinematic condition,
$$\begin{eqnarray}v_{i}=S_{t}+u_{i}S_{x},\quad i=1,2,\end{eqnarray}$$
and continuity of normal and tangential stresses,
$$\begin{eqnarray}\displaystyle & \displaystyle \left[-p_{i}(1+S_{x}^{2})+\frac{2}{\mathit{Re}_{i}}(S_{x}^{2}u_{ix}+v_{iy}-S_{x}(u_{iy}+v_{ix}))\right]_{2}^{1}=\frac{1}{\mathit{We}}\frac{S_{xx}}{\sqrt{1+S_{x}^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[4\frac{\unicode[STIX]{x1D707}_{i}}{\unicode[STIX]{x1D707}_{1}}\,S_{x}u_{ix}+\frac{\unicode[STIX]{x1D707}_{i}}{\unicode[STIX]{x1D707}_{1}}(S_{x}^{2}-1)(u_{iy}+v_{ix})\right]_{2}^{1}=0, & \displaystyle\end{eqnarray}$$
$[f_{i}]_{2}^{1}=f_{1}-f_{2}$
is used. Here,
$\mathit{We}=\unicode[STIX]{x1D70C}U^{2}d/\unicode[STIX]{x1D6FE}$
is the Weber number, also equal to
$\mathit{We}=\mathit{Ca}\mathit{Re}_{2}$
, with
$\mathit{Ca}=\unicode[STIX]{x1D707}_{2}U/\unicode[STIX]{x1D6FE}$
the capillary number.
When the interface is flat,
$y=\ell$
, the system (3.1)–(3.6) admits an exact base state solution
$\boldsymbol{u}_{1}=(\overline{U}_{1}(y),0)$
,
$\boldsymbol{u}_{2}=(\overline{U}_{2}(y),0)$
, where
$$\begin{eqnarray}\overline{U}_{1}(y)=\frac{my}{\ell (m-1)+1},\quad \overline{U}_{2}(y)=1-\frac{1-y}{\ell (m-1)+1}.\end{eqnarray}$$
One of the main objectives of the present work is to study the nonlinear stability of this state, and in what follows analytical progress is made for thin upper layers – the theory is also evaluated against DNS and experiments.
Taking the mean thickness of region 2 to be
$\unicode[STIX]{x1D716}\ll 1$
(analogously
$\ell =1-\unicode[STIX]{x1D716}$
), we introduce weakly nonlinear perturbations
$S(x,t)=(1-\unicode[STIX]{x1D716})-\unicode[STIX]{x1D716}^{2}\widetilde{H}(x,t)$
, with
$\widetilde{H}=O(1)$
. The canonical scaling
$\mathit{Ca}=\unicode[STIX]{x1D716}\mathit{Ca}_{0}$
, with
$\mathit{Ca}_{0}=O(1)$
, is introduced to ensure coupling between the two layers. The analysis follows the work of Kalogirou & Papageorgiou (Reference Kalogirou and Papageorgiou2016) (see also Kalogirou Reference Kalogirou2014), the main difference being that in the present case the thin layer is in the vicinity of the upper moving plate to enable direct comparisons with experiments. The expansions in region 1 are
$u_{1}=\overline{U}_{1}(y)+\unicode[STIX]{x1D716}^{2}U_{1}(x,y,t)+\cdots$
,
$v_{1}=\unicode[STIX]{x1D716}^{2}V_{1}(x,y,t)+\cdots$
,
$p_{1}=\unicode[STIX]{x1D716}^{2}P_{1}(x,y,t)+\cdots$
; in region 2 we introduce a stretched variable
$\unicode[STIX]{x1D701}=(1-y)/\unicode[STIX]{x1D716}$
(the undisturbed interface is now at
$\unicode[STIX]{x1D701}=1$
) and write
$u_{2}=\overline{U}_{2}(\unicode[STIX]{x1D701};\unicode[STIX]{x1D716})+\unicode[STIX]{x1D716}^{3}U_{2}(x,\unicode[STIX]{x1D701},t)+\cdots$
,
$v_{2}=\unicode[STIX]{x1D716}^{4}V_{2}(x,\unicode[STIX]{x1D701},t)+\cdots$
,
$p_{2}=\unicode[STIX]{x1D716}P_{2}(x,\unicode[STIX]{x1D701},t)+\cdots$
. The solutions then retain a leading-order balance between the pressure gradient and viscous terms in the momentum equations.
At leading order we have advection with the undisturbed interfacial velocity
$U_{s}=U_{2}|_{y=1-\unicode[STIX]{x1D716}}=U_{1}|_{\unicode[STIX]{x1D701}=1}$
, as expected. An evolution equation for
$\widetilde{H}$
is found at the next order after transforming to a Galilean frame of reference and introducing a slow time scale
$$\begin{eqnarray}\displaystyle & \displaystyle X=x-U_{s}\,t,\quad \tilde{t}=\unicode[STIX]{x1D716}^{2}\,t, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left.\widetilde{H}_{\tilde{t}}-\frac{1}{m}\widetilde{H}\widetilde{H}_{X}+\frac{1}{3\mathit{Ca}_{0}}\widetilde{H}_{XXXX}-\frac{1}{2m}{\mathcal{T}}\right|_{y=1}=0. & \displaystyle\end{eqnarray}$$
In addition to unsteadiness and nonlinearity, equation (3.9) contains surface tension and involves the non-local inertial term
${\mathcal{T}}(x,y)=U_{1xy}+V_{1xx}$
which is found by solving in region 1, as explained next. The slow time dynamics and the
$O(\unicode[STIX]{x1D716}^{2})$
perturbations imply that the flow in region 1 is governed by the linearised steady Navier–Stokes equations at
$\mathit{Re}_{1}=O(1)$
. Taking Fourier transforms, eliminating the pressure
$P_{2}$
and writing
$\widehat{V}_{1}=-\text{i}k(1/m-1)\widehat{\widetilde{H}}F(y)$
(hats denote Fourier transforms), yields the following Orr–Sommerfeld-type problem:
$$\begin{eqnarray}\displaystyle & (F^{\prime \prime \prime \prime }-2k^{2}F^{\prime \prime }+k^{4}F)-\text{i}k\,\mathit{Re}_{1}\,y(F^{\prime \prime }-k^{2}F)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & F(0)=0,\quad F^{\prime }(0)=0,\quad F(1)=0,\quad F^{\prime }(1)=1. & \displaystyle\end{eqnarray}$$
$y=0$
and continuity of velocities at the interface
$y=1$
. It follows that the non-local term in (3.9) is given by 
$$\begin{eqnarray}{\mathcal{T}}|_{y=1}=\frac{\text{i}}{\unicode[STIX]{x03C0}}\left(1-\frac{1}{m}\right)\int _{-\infty }^{+\infty }{\mathcal{N}}(k;\mathit{Re}_{1}){\hat{H}}(k)\text{e}^{\text{i}kx}\,\text{d}k,\quad {\mathcal{N}}(k;\mathit{Re}_{1})=-\frac{k}{2}F^{\prime \prime }(1).\end{eqnarray}$$
The final rescalings
$\widetilde{H}=(m/3Ca_{0})H$
,
$\widetilde{t}=(3Ca_{0}/m)T$
produce the canonical equation
$$\begin{eqnarray}H_{T}-HH_{X}+H_{XXXX}-\frac{\text{i}\unicode[STIX]{x1D6EC}}{2\unicode[STIX]{x03C0}}\int _{-\infty }^{+\infty }{\mathcal{N}}(k;\mathit{Re}_{1})\widehat{H}(k)\text{e}^{\text{i}kX}\,\text{d}k=0.\end{eqnarray}$$
The parameter
$\unicode[STIX]{x1D6EC}=3\mathit{Ca}_{0}(1/m)(1-1/m)$
represents the effects of viscosity stratification. Linearising about
$H=0$
and looking for solutions proportional to
$\text{e}^{\text{i}kx+\unicode[STIX]{x1D70E}t}$
yields
$$\begin{eqnarray}\unicode[STIX]{x1D70E}=-k^{4}-\frac{\text{i}k\unicode[STIX]{x1D6EC}}{4\unicode[STIX]{x03C0}}F^{\prime \prime }(1).\end{eqnarray}$$
Instability is found if
$\text{Re}(\unicode[STIX]{x1D70E})=(k/4\unicode[STIX]{x03C0})\unicode[STIX]{x1D6EC}\,\text{Im}(F^{\prime \prime }(1))-k^{4}>0$
, where
$\text{Re}$
and
$\text{Im}$
denote the real and imaginary parts of a given quantity. Numerical solution of the boundary value problem (3.10) predicts that
$\text{Im}(F^{\prime \prime }(1))\geqslant 0$
(we do not have a proof for this). Therefore, the flow is unstable only if the film is more viscous than the lower fluid (
$\unicode[STIX]{x1D6EC}>0$
), completely in line with the results of Yih (Reference Yih1967).
Equation (3.12) is solved numerically on
$2L$
-periodic domains. The non-local term becomes a Fourier series, and the initial value problem to solve is
$$\begin{eqnarray}\displaystyle & \displaystyle H_{T}-HH_{X}+H_{XXXX}-\text{i}\unicode[STIX]{x1D6EC}\mathop{\sum }_{k=-\infty }^{+\infty }{\mathcal{N}}\left(\frac{k\unicode[STIX]{x03C0}}{L};\mathit{Re}_{1}\right)\widehat{H}(k)\text{e}^{\text{i}k\unicode[STIX]{x03C0}X/L}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle H(X,0)=A\sin \left(\frac{\unicode[STIX]{x03C0}X}{L}\right), & \displaystyle\end{eqnarray}$$
where
$A>0$
is an amplitude set to 1 unless stated otherwise. Other initial conditions can be used – for example
$H(X,0)=A\sin (2\unicode[STIX]{x03C0}X/L)$
when searching for bistability. Equation (3.14) is discretised spectrally in space, and implicit–explicit BDF (backward difference) schemes are used for the time-stepping. The numerical results that follow were obtained via a second-order-accurate time-discretisation scheme; for numerical analysis details including convergence theorems about BDF schemes (with as high as sixth-order accuracy) see Akrivis, Papageorgiou & Smyrlis (Reference Akrivis, Papageorgiou and Smyrlis2012).
We conclude this section by constructing the data that are used in experimental and DNS comparisons discussed in detail in § 4. The experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995) are in the inertial regime (
$Re_{1}\approx 500$
), and identify, among other phenomena, interfacial travelling waves of permanent form. The interfacial amplitude is recorded by a fixed probe, and the spatial profiles follow from this signal if the wave speed is known – this is also implemented in computations. The value of the spatial period
$2L$
is set by the experiment (from table 1 we have
$2L=20\unicode[STIX]{x03C0}\approx 63$
), and at such values (3.14) has travelling wave solutions of the form
$H(X,T)\equiv H(X-cT)$
, where
$c$
is the constant phase velocity. In model computations, we record the amplitude at
$X=0$
, and so the period of oscillation in time is given by
$2L/c$
– this needs to be multiplied by
$c/U_{s}$
in order to compare with experiments and DNS which place a probe at
$x=0$
and not
$x=U_{s}t$
. In a direct numerical simulation for a given finite
$\unicode[STIX]{x1D716}$
, we monitor
$$\begin{eqnarray}\displaystyle S(0,t)=1-\unicode[STIX]{x1D716}-f(t), & & \displaystyle\end{eqnarray}$$
where
$f(t)$
is periodic (for small
$\unicode[STIX]{x1D716}$
, the period is approximately
$2L/U_{s}$
). We scale the oscillatory part as in the experiments by defining
$S(t)=S(0,t)-(1-\unicode[STIX]{x1D716})$
and use the numerical data to find
$A_{+}=\max \{S(t)\}$
and
$A_{-}=\min \{S(t)\}$
. The scaled amplitude is
$$\begin{eqnarray}\displaystyle \frac{S(t)}{S_{sat}},\quad \text{where}~S_{sat}=\frac{1}{2}(A_{+}+A_{-}). & & \displaystyle\end{eqnarray}$$
An analogous amplitude normalisation is used in computations of the model also, before experimental or DNS comparisons are carried out.
4 Comparison between DNS, model computations and experiments
In this section, we compare computed solutions of (3.14) with those coming from DNS, and use the model solutions to make quantitative comparisons with experimental results (the experiments are on very long domains which are challenging to simulate directly).
4.1 Direct numerical simulations and comparison with model solutions
We carry out DNS of the flow configuration given in figure 1 and the fluids described in table 1 as reported in the experiments. The initial condition in each fluid layer is a superposition of the basic states (3.7) and a perturbation (chosen to match initial condition (3.15)) of the form
$$\begin{eqnarray}S(x,0)=1-\unicode[STIX]{x1D716}-\unicode[STIX]{x1D716}^{2}\sin \left(\frac{\unicode[STIX]{x03C0}x}{L}\right),\end{eqnarray}$$
with domain lengths
$2L=4,8,12$
.
The problem was implemented using the volume-of-fluid (VOF) open-source software Gerris (Popinet Reference Popinet2003). The quadtree-based structure of the code enabled adaptive mesh refinement and efficient parallelisation; the numerical schemes employed are second-order accurate in both space and time. Periodic boxes representing the three geometries are constructed in Gerris, and each sub-unit square is discretised using a uniform grid with 64 nodes in each dimension; two additional refinement levels are prescribed around the fluid–fluid interface so that the local grid spacing here becomes
$1/4$
of the bulk mesh. For the largest and most expensive geometries, a total of at least 12 288 cells were used to represent the interfacial vicinity, while the total number of degrees of freedom in the domain was approximately
$6\times 10^{4}$
. All simulations were run for sufficiently large times to enable the interfacial profiles to reach saturation to travelling waves. In each case, this required several days of runtime on local high-performance computing facilities, with calculations executed in parallel on 4, 8 and 12 CPUs respectively.
Numerical computations and comparisons between DNS and the model are presented in figure 2 for
$\mathit{Re}_{1}=516$
, corresponding to experiment C in table 2 (the other parameters
$\mathit{Re}_{2}$
,
$We$
and
$Ca$
are as given in the table). The interfacial amplitude is monitored at
$x=0$
, the centre of the domain, and normalised as described in § 3 – see (3.16) and (3.17). Computations of the evolution equation (3.14)–(3.15) were also carried out, with the value of
$\unicode[STIX]{x1D6EC}$
adjusted in order to obtain the best possible agreement; it should be noted that the only adjustable parameter in
$\unicode[STIX]{x1D6EC}$
is
$\mathit{Ca}_{0}$
, where
$\mathit{Ca}=\unicode[STIX]{x1D716}\,\mathit{Ca}_{0}$
.
Figure 2. (a,c) Comparison of wave periods (in seconds) computed with the model (diamonds) and DNS (circles) as a function of
$\unicode[STIX]{x1D716}$
varying from 0.01 to 0.2. The Reynolds number is
$\mathit{Re}_{1}=516$
, and two channels with different aspect ratios are used (4 : 1 (a), 8 : 1 (c)). The period obtained by tracking a particle on the undisturbed interface is included for reference (filled squares). (b,d) Position of the interface at
$x=0$
obtained numerically using the model (solid curves) and DNS (dashed lines). (b) Aspect ratio 4 : 1 and
$\unicode[STIX]{x1D716}=0.1$
; (d) aspect ratio 8 : 1 and
$\unicode[STIX]{x1D716}=0.2$
.
A systematic comparison between model and DNS is given in figure 2 for channel aspect ratios of 4 : 1 and 8 : 1. Panels (a,c) show the computed periods of oscillation in seconds versus the upper layer thickness
$\unicode[STIX]{x1D716}$
– open circles are used for the DNS, diamonds for the model and filled squares represent the time period of a particle on the undisturbed interface with velocity given by (3.7). The wave speed is found to be independent of the channel length, and we have confirmed that as
$\unicode[STIX]{x1D716}$
becomes smaller it tends to the undisturbed base velocity (this is asymptotically equal to the upper plate speed
$U=0.244~\text{m}~\text{s}^{-1}$
). We also find that the period of oscillation increases proportionally to the channel length – for the
$8\times 1$
channel the period is approximately double that of the
$4\times 1$
channel (to within an error of
$O(\unicode[STIX]{x1D716}^{2})$
). We note that the periods predicted by the model slightly underestimate the values obtained by DNS.
A direct comparison of the travelling wave shapes is given in figure 2(b,d) for the 4 : 1 and 8 : 1 aspect ratio geometries respectively. The DNS results are depicted with a dashed red curve and the model ones with a solid curve – the agreement is seen to be very good. The model computations were obtained by setting
$\unicode[STIX]{x1D6EC}=0.4$
, implying that
$\mathit{Ca}_{0}\approx 0.577$
(the value of
$m=2.76$
is fixed by the experiments). The simulations take
$\mathit{Ca}=0.242$
, and using
$\unicode[STIX]{x1D716}=0.1$
predicts
$\mathit{Ca}_{0}\approx 2.42$
, which is larger than the value used in the model. This discrepancy is due to the asymptotic nature of the evolution equations and we do not expect an exact correspondence at finite
$\unicode[STIX]{x1D716}$
. In the longer domain of figure 2(d), the dimensionless film thickness is set to
$\unicode[STIX]{x1D716}=0.2$
and the model uses
$\unicode[STIX]{x1D6EC}=0.2$
and
$\mathit{Ca}_{0}=0.29$
. Once again, the agreement between DNS results (red dashed curve) and model predictions (black solid curve) is very good. More notably, the interfacial shapes in the model–DNS comparisons have the same characteristics, with features due to the nonlinearity dominating as the domain size increases. Finally, we note that for both channel aspect ratios, the cases presented correspond to the largest values of
$\unicode[STIX]{x1D716}$
, thus providing the most challenging conditions for the comparison.
4.2 Comparisons with the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995)
Given the success of the model in comparisons with DNS, we turn next to comparisons with the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995). As noted earlier, these require relatively long domains (approximately 63 dimensionless units), which are challenging for DNS. In what follows, we compare with experiments A, B and C in table 2, corresponding to
$\mathit{Re}_{1}=389$
, 461 and 516 respectively. For the model computations, we set
$2L=20\unicode[STIX]{x03C0}\approx 62.832$
, and calculate the required kernels (3.11) at the different values of
$\mathit{Re}_{1}$
. The parameter
$\unicode[STIX]{x1D6EC}$
is set to
$\unicode[STIX]{x1D6EC}=0.0015$
in order to retain a sufficient number of unstable modes to enable dynamical comparisons. Signals are constructed as described at the end of § 3 and the results are presented in figure 3, which reproduces the experimental results in (a) and the computational ones in (b). The wave amplitudes are normalised with the saturated amplitude of the wave for the largest Reynolds number
$\mathit{Re}_{1}=516$
, as done in the experiments. In (b), we also superimpose with dotted curves numerical results obtained by solving the localised equation
$$\begin{eqnarray}H_{T}-HH_{X}+\unicode[STIX]{x1D6FC}H_{XXXX}+2\unicode[STIX]{x1D6EC}H_{X}+\frac{\unicode[STIX]{x1D6EC}\,\mathit{Re}_{1}}{20}H_{XX}-\unicode[STIX]{x1D6EC}\left(\frac{2}{15}+\frac{17\,\mathit{Re}_{1}}{50\,400}\right)H_{XXX}=0,\end{eqnarray}$$
which arises by taking the long-wave limit in (3.14) and introducing a ‘correction’ parameter
$\unicode[STIX]{x1D6FC}=1.125$
so that the non-local and local models have the same band of linearly unstable modes (with this choice, the maximum growth rate of the local model is 0.01195, compared with 0.01170 for that of the non-local model). The numerical computations indicate that the non-local model does better in comparisons with the experiments. The reason for including results from (4.2) is to motivate a possible rational way of constructing appropriate dispersive Kuramoto–Sivashinsky-like equations whose dynamics resemble experimental observations – e.g. see Barthelet & Charru (Reference Barthelet, Charru, Renardy, Coward, Papageorgiou and Sun1995). The use of a local model is not commendable since it nevertheless requires the analysis of the non-local equation (3.14). The latter is straightforward to implement numerically and hence preferable. The evolution of the amplitudes of the first three harmonics of the non-local and local models was also compared for the parameters of figure 3. We find (results not shown) that for
$\mathit{Re}_{1}=389$
the local model underpredicts the amplitudes, whereas for the largest value
$\mathit{Re}_{1}=516$
it overpredicts them.
Figure 3. Position of the interface at a fixed point. (a) Experimental shapes obtained by Barthelet et al. (taken from Barthelet et al. (Reference Barthelet, Charru and Fabre1995), p. 36) and (b) interfacial shapes obtained numerically by solving the model equation (3.14) and its localised version (4.2). The three plots correspond to the values of upper plate speeds or Reynolds numbers given in table 2.
The waves in figures 3(a) and 3(b) have similar shapes, consisting of steep wavefronts and sharp troughs which are enhanced as
$\mathit{Re}_{1}$
increases. The period of oscillation decreases as
$\mathit{Re}_{1}$
increases, and a comparison between theory and experiment is provided in table 3. The difference ranges from approximately 10 % for experiment A to 14 % for experiment C. The model underestimates the experimental values, and this is consistent with the fact that the experiments study stably stratified flows (the density ratio is approximately 0.74 – see table 1) rather than equal densities. In fact, preliminary analysis indicates that stable density stratification modifies (3.14) with a second-order diffusion term that would increase the period of oscillations, all other parameters being equal.
Table 3. Comparison between periods of oscillation for experiments A, B, C given in table 2 and the corresponding ones from the model computations. All values are given in seconds.
4.3 Predicting bistability using the model
The experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995) find bistable travelling waves (of wavelengths equal to the channel perimeter or half the perimeter respectively) for sufficiently large values of the upper plate speed – see their figures 22–24. This phenomenon is reproduced by the model equation, as we describe next. In particular, we model the experiment having upper plate velocity
$U=0.335~\text{m}~\text{s}^{-1}$
, implying a Reynolds number of
$\mathit{Re}_{1}\approx 709$
. We use
$\unicode[STIX]{x1D6EC}=0.002$
so that two unstable modes exist, and take initial conditions of the form of either
$H(X,0)=0.005\sin (\unicode[STIX]{x03C0}X/L)$
or
$H(X,0)=0.005\sin (2\unicode[STIX]{x03C0}X/L)$
respectively to excite these modes independently to nonlinear bistable saturated states. The results are shown in figure 4, with the saturated profiles on the left and the evolution of the corresponding first four harmonics on the right. Saturation takes place fairly rapidly for the
$2L$
-periodic waves, and the saturated amplitudes of the harmonics decrease as
$k$
increases. For the
$L$
-periodic wave (c,d), saturation takes longer (approximately
$6T_{sat}$
s, with
$T_{sat}=4.68$
s being the oscillation period after saturation). In addition, the bimodal nature is preserved throughout the evolution, with only the even modes
$k=2$
and
$k=4$
being present. It is interesting to note that the
$L$
-periodic solutions (c) are almost symmetric, in contrast to the more pronounced asymmetry observed for
$2L$
-periodic ones (a) and the experiments of Barthelet et al. (Reference Barthelet, Charru and Fabre1995) (see figure 22 of that paper). A more detailed study of the dynamical system which fully explores the solution phase space and tracks bifurcations is warranted and is left for future work.
Figure 4. Evolution of the interfacial position at a fixed point (a,c), and amplitude of harmonics (b,d) showing bistability. Parameter values:
$\mathit{Re}_{1}=709$
,
$\unicode[STIX]{x1D6EC}=0.002$
.
$T_{sat}=4.68$
s is the saturated wave period. The prescribed initial conditions are
$H(X,0)=0.005\sin (\unicode[STIX]{x03C0}X/L)$
(a,b) and
$H(X,0)=0.005\sin (2\unicode[STIX]{x03C0}X/L)$
(c,d).
5 Discussion
A novel evolution equation has been proposed and studied to describe the nonlinear stability of two-layer Couette flows. The equation is first validated using DNS of the Navier–Stokes equations, and then used to make comparisons with appropriate experiments from the work of Barthelet et al. (Reference Barthelet, Charru and Fabre1995). The analysis is performed for the thickness
$\unicode[STIX]{x1D716}$
of the layer adjacent to the moving plate being small. The smallest thicknesses available in the experiments are
$\unicode[STIX]{x1D716}=0.2$
, and, even though this is not necessarily small, the model does fairly well in reproducing experimental observations at Reynolds numbers as large as
$O(10^{3})$
. A crucial aspect of our model is that it fully captures inertia in the thicker layer through a linear non-local term whose spectral properties are found by solving an Orr–Sommerfeld-type problem. As a result, we can set the Reynolds number
$\mathit{Re}_{1}$
and viscosity ratio
$m$
directly from an experiment, and as shown in § 3 this leaves us with the parameter
$\unicode[STIX]{x1D6EC}=3\mathit{Ca}_{0}(1/m)(1-1/m)$
, where
$Ca=\unicode[STIX]{x1D707}_{2}U/\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D716}\,Ca_{0}$
is the canonical scaling for the capillary number which retains coupling between the two fluid regions. In a given experiment or direct numerical computation, one can calculate the value of
$Ca$
– see table 2 for example – to ensure that it is small. The value of
$Ca_{0}$
cannot be found from such data, and the parameter
$\unicode[STIX]{x1D6EC}$
is not fixed. Since the wave periods are fixed by the experiment,
$\unicode[STIX]{x1D6EC}$
is the only parameter that must be selected in making comparisons with experiments and DNS (note that the functional form of the linear spectrum is not affected by
$\unicode[STIX]{x1D6EC}$
since it enters as a multiplicative factor – see (3.13)). The model is computed and
$\unicode[STIX]{x1D6EC}$
is modified until agreement is obtained, as shown in the results. If
$\unicode[STIX]{x1D6EC}$
is too small or too large, the dynamics can range from trivial to chaotic, and hence is at variance with DNS and experiments.
The model has also been shown to predict bistability, with unimodal and bimodal solutions coexisting at the same parameter values. This has been demonstrated for the experimental parameters of figure 22 of Barthelet et al. (Reference Barthelet, Charru and Fabre1995), and the results are described in § 4.3. Other phenomenologically motivated model equations – see Barthelet & Charru (Reference Barthelet, Charru, Renardy, Coward, Papageorgiou and Sun1995) – have been shown to predict bistable steady states, but it is difficult to make comparisons with experiments, unlike for the inertia-retaining models derived and studied here.
As far as we know, this work is the first to compare directly solutions of model equations with those of DNS. The latter take orders of magnitude more time to run even on state-of-the-art high-performance computers. The excellent agreement demonstrated makes such models good candidates in the description of dynamics beyond the capabilities of standard lubrication-type approximations of multi-fluid flows – see Craster & Matar (Reference Craster and Matar2009). Our novel methodology can also be extended to viscoelastic flows in order to assess the effects of elasticity on the stability and nonlinear dynamics arising from the interaction between thick and thin layers. Such extensions are the subject of ongoing work.
Acknowledgements
D.T.P. and R.C. were partly supported by EPSRC grants EP/K041134 and EP/L020564. A.K. acknowledges the support of a Roth PhD scholarship by the Department of Mathematics, Imperial College London.
