1 Introduction
The displacement of one fluid by another in confined geometry is one of the most common processes in nature. These interesting flows also find numerous industrial applications in the petroleum industry (oil and gas exploration, extraction and transportation; see e.g. Nelson & Guillot Reference Nelson and Guillot2006), in manufacturing (coating flows, co-extrusion; see e.g. Schweizer & Kistler Reference Schweizer and Kistler2012), in biomedical applications (digestion, clearing of mucus plugs in alveoli; see e.g. Huh et al. Reference Huh, Fujioka, Tung, Futai, Paine, Grotberg and Takayama2007), and in food processing (biofilm removal, cleaning of processing machinery; see e.g. Wiklund, Stading & Trägårdh Reference Wiklund, Stading and Trägårdh2010). Displacement flows may emerge in uniform geometries but their occurrence in geometries with at least a degree of non-uniformity is much more prevalent (Al-Housseiny, Tsai & Stone Reference Al-Housseiny, Tsai and Stone2012). Through developing a semi-analytical model, our aim in this paper is to provide a basic understanding of the effects of the simplest geometrical non-uniformity on laminar displacement flows.
In this work, we consider miscible displacement flows along a slightly non-uniform two-dimensional plane channel of width
$\hat{D}_{0}+\tan {\it\alpha}\hat{x}$
(i.e.
${\approx}\hat{D}_{0}+{\it\alpha}\hat{x}$
considering
$|{\it\alpha}|\ll 1$
; cf. figure 1) in the large-Péclet-number limit,
$Pe=\hat{D}_{0}\hat{V}_{0}/\hat{D}_{m}\gg 1$
, where
$\hat{D}_{m}$
is the molecular diffusivity and
$\hat{V}_{0}$
is the mean displacement velocity. The consideration of this limit implies that the two fluids do not have sufficient time to mix on the time scales of interest. In this case, instead of a concentration diffusion-equation approach, we rely on a kinematic equation as a reasonable approximation to model the two fluids. This is in accordance with the previous studies finding that the high-Péclet-number regime approaches the zero-surface-tension immiscible limit (
$Pe\rightarrow \infty$
), provided that the flow remains stable (see e.g. Chen & Meiburg Reference Chen and Meiburg1996; Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Rakotomalala, Salin & Watzky Reference Rakotomalala, Salin and Watzky1997; Yang & Yortsos Reference Yang and Yortsos1997). In the limit of large-
$Pe$
displacement flows, we consider the flows with the following characteristics: (1) buoyancy is a significant driving force affecting the fluids’ motion; (2) the channel is slightly converging or diverging by simply assuming that the channel upper wall is slightly inclined with respect to the lower wall; (3) a viscous flow regime can be found for at least relatively small Reynolds numbers (defined later); and (4) the lubrication/thin-film approximation is the principal tool in this study and the scenarios that we consider involve fluids with Newtonian or non-Newtonian Herschel–Bulkley rheology.
Miscible displacement flows of Newtonian fluids in uniform pipes and channels have been studied in detail. For example, once viscous effects are negligible, Benjamin (Reference Benjamin1968), Shin, Dalziel & Linden (Reference Shin, Dalziel and Linden2004) and Birman et al. (Reference Birman, Battandier, Meiburg and Linden2007) studied displacement flows that are characterized by a balance between buoyancy and internal stresses. Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2010, Reference Taghavi, Seon, Wielage-Burchard, Martinez and Frigaard2011, Reference Taghavi, Alba, Seon, Wielage-Burchard, Martinez and Frigaard2012c ), Taghavi, Alba & Frigaard (Reference Taghavi, Alba and Frigaard2012a ) and Alba, Taghavi & Frigaard (Reference Alba, Taghavi and Frigaard2012) studied density-unstable and density-stable buoyant miscible displacement flows along inclined pipes and channels and addressed the effects of pipe/channel inclinations, viscosity and density ratios and the imposed flow rate. Seon et al. (Reference Seon, Hulin, Salin, Perrin and Hinch2004, Reference Seon, Hulin, Salin, Perrin and Hinch2005, Reference Seon, Hulin, Salin, Perrin and Hinch2006, Reference Seon, Znaien, Salin, Hulin, Hinch and Perrin2007b ) studied in depth the exchange flow counterparts, i.e. with a zero imposed flow. Many authors have also considered Newtonian displacement flows in capillary tubes (Chen & Meiburg Reference Chen and Meiburg1996; Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Vanaparthya & Meiburg Reference Vanaparthya and Meiburg2008) as well as in horizontal and vertical Hele-Shaw (HS) geometries (see e.g. Talon, Goyal & Meiburg Reference Talon, Goyal and Meiburg2013; Heussler et al. Reference Heussler, Oliveira, John and Meiburg2014, respectively).
Displacement flows of non-Newtonian fluids have also received attention. However, this has been mostly allied to studying the well-known viscous fingering (VF) instability and pattern formations in HS geometries, for which of course significant developments date back to the classical stability analysis of Saffman & Taylor (Reference Saffman and Taylor1958) considering the modification of the gap-averaged HS theory where there is a wetting film of the displaced fluid at the walls. A whole diverse class of problems has been studied in this domain for non-Newtonian fluids; see e.g. various features studied in Coussot (Reference Coussot1999), Fast et al. (Reference Fast, Kondic, Shelley and Palffy-Muhoray2001), Nguyen et al. (Reference Nguyen, Folch, Verma, Henry and Plapp2010) and Ebrahimi, Taghavi & Sadeghy (Reference Ebrahimi, Taghavi and Sadeghy2015), snowflake-like patterns in Buka, Palffy-Muhoray & Racz (Reference Buka, Palffy-Muhoray and Racz1987), and branched, fractal, or fracture-like structures in Nittmann, Daccord & Stanley (Reference Nittmann, Daccord and Stanley1985), Lemaire et al. (Reference Lemaire, Levitz, Daccord and Van Damme1991) and Ignes-Mullol, Zhao & Maher (Reference Ignes-Mullol, Zhao and Maher1995). Shear-thinning properties may induce side branching or crack-like patterns (Kondic, Shelley & Palffy-Muhoray Reference Kondic, Shelley and Palffy-Muhoray1998; Ben Amar & Poire Reference Ben Amar and Poire1999) and yield-stress features strikingly modify the morphological patterns (Lindner, Coussot & Bonn Reference Lindner, Coussot and Bonn2000).
The literature on non-Newtonian displacement flows in geometries other than HS is less developed. Bittleston, Ferguson & Frigaard (Reference Bittleston, Ferguson and Frigaard2002) derived a two-dimensional model of laminar non-Newtonian displacement flows in an annulus, with applications in the primary cementing of oil and gas wells, and they showed that a channel of viscoplastic mud may be left behind in the narrow side of the annulus. Dimakopoulos & Tsamopoulos (Reference Dimakopoulos and Tsamopoulos2003, Reference Dimakopoulos and Tsamopoulos2007) and De Sousa et al. (Reference De Sousa, Soares, de Queiroz and Thompson2007) investigated gas–liquid displacements in tubes, finding residual layers in steady-state displacements. Zhang & Frigaard (Reference Zhang and Frigaard2006) studied displacement flow for a range of simple non-Newtonian fluids. Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009) studied heavy–light and light–heavy displacement of viscoplastic fluids in a channel close to horizontal. They showed that a yield stress in the displacing fluid enhances the displacement flow process, while a yield stress in the displaced fluid leads to completely static residual wall layers. These interesting layers have been first analysed in the displacement flow context by Allouche, Frigaard & Sona (Reference Allouche, Frigaard and Sona2000). By considering a symmetric displacement flow of two Bingham fluids in a plane channel, they showed that, for these layers to exist, the yield stress of the displaced fluid has to exceed that of the displacing fluid. Their work was extended by Wielage-Burchard & Frigaard (Reference Wielage-Burchard and Frigaard2011), whose computational study provided further insight into the effects of the dimensionless parameters (i.e. the Reynolds number, Bingham number and viscosity ratio) and included the effects of flow rate oscillations. Taghavi et al. (Reference Taghavi, Alba, Moyers-Gonzalez and Frigaard2012b ) and Alba et al. (Reference Alba, Taghavi, de Bruyn and Frigaard2013a ) presented experimental results of a viscoplastic fluid displacement flow in inclined pipes. They focused on the situations where the yield stress is much larger than a typical viscous stress in the displacing fluid and identified a ‘central-type’ and a ‘slump-type’ displacement regime for higher density ratios.
Fluid flows along non-uniform channels have been studied over many years (e.g. Hasegawa & Izuchi Reference Hasegawa and Izuchi1983; Nishimura et al. Reference Nishimura, Arakawa, Murakami and Kawamura1989), due to their industrial applications, e.g. increasing heat transfer (Oviedo-Tolentino et al. Reference Oviedo-Tolentino, Romero-Méndez, Hernández-Guerrero and Girón-Palomares2008), peristaltic pumping (Burns & Parkes Reference Burns and Parkes1967), etc. Displacements fluid flows in non-uniform channels have also recently received attention. Recent studies (e.g. Al-Housseiny et al. Reference Al-Housseiny, Tsai and Stone2012; Wilson Reference Wilson2012; Al-Housseiny & Stone Reference Al-Housseiny and Stone2013; Dias & Miranda Reference Dias and Miranda2013; Al-Housseiny Reference Al-Housseiny2014), have considered immiscible Newtonian displacements in a tapered rectilinear HS by introducing a small depth gradient, crafting a slightly converging or diverging HS and showing that certain flow instabilities (e.g. VF) can be effectively stabilized through a variable transverse curvature mechanism. Other relevant works (e.g. Pihler-Puzović et al. Reference Pihler-Puzovic, Illien, Heil and Juel2012, Reference Pihler-Puzović, Périllat, Russell, Juel and Heil2013; Al-Housseiny, Christov & Stone Reference Al-Housseiny, Christov and Stone2013) have concentrated on Newtonian displacements in a radial HS, finding that VF can be fascinatingly suppressed when the upper boundary of the cell was replaced by an elastic membrane. They have observed that the elastic wall creates a non-uniform displacement flow geometry that strongly influences the finger shape and its development, replacing the characteristic dendritic patterns in a traditional HS with a large number of very short fingers.
Relevant to our topic is a well-known argument, referred to as the ‘lubrication paradox’, first made by Lipscomb & Denn (Reference Lipscomb and Denn1984), who argued that the use of classical lubrication scaling techniques for yield-stress fluids in spatially varying geometries results in finding plug speed that slowly varies in the main flow direction. This slow variation of the plug speed implies that plug regions predicted by classical lubrication models are slowly deformable and therefore they are not ‘truly’ rigid. This is in contrast to analytical and computational studies, which have confirmed the existence of true unyielded plug regions (which are not deformable) in complex geometries (Putz, Frigaard & Martinez Reference Putz, Frigaard and Martinez2009). Walton & Bittleston (Reference Walton and Bittleston1991) and Putz et al. (Reference Putz, Frigaard and Martinez2009) defined pseudo-plug regions while considering the problematic deformable plugs for which the shear stress slightly exceeds the yield stress. The observation of the pseudo-plugs has driven new discussions in the literature, e.g. on the possibility of truly rigid plugs within the pseudo-plugs or the development of pseudo-plugs to true plugs for sufficiently weak spatial variations (see Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014, for a review). Frigaard & Ryan (Reference Frigaard and Ryan2004) developed an asymptotic solution for the Poiseuille flow of a Bingham fluid along a channel of slowly varying width with a small-amplitude long-wavelength perturbation (
$a$
). They showed that an asymptotic solution with an intact unyielded plug region can be found for
$a<a_{c}$
(critical) with
$a_{c}\sim O({\it\delta})$
, where
${\it\delta}$
is the aspect ratio. Lately, Roustaei & Frigaard (Reference Roustaei and Frigaard2013) extended the work mentioned to larger amplitude perturbations. Relevant to our work, we emphasize that there is no lubrication paradox. In fact, identifying where there are true unyielded regions rather than pseudo-unyielded regions for displacement of two viscoplastic fluids only requires a more elaborate asymptotic analysis (Balmforth et al.
Reference Balmforth, Frigaard and Ovarlez2014). In the current paper, our goal is not to resolve these details by developing a higher-order lubrication model; instead, we would like to provide a basic understanding of displacement flows in non-uniform geometry based on the leading-order lubrication theory. Therefore, the predictions of our model for the cases with a yield-stress rheology, focused on the appearance of static residual wall layers (see § 5.2.1), may be qualitative.
As discussed, the literature addressing non-Newtonian displacement flows in non-uniform geometries is not well developed, owing to the challenges in this area. Firstly, these flows are very complex to study in generality because of the number and the range of the dimensionless parameters involved. This has limited the literature offering comprehensive pictures of various flow regimes and predictions of the main flow features. In addition, despite some efforts (e.g. Hallez & Magnaudet Reference Hallez and Magnaudet2008), the role of geometry in general and the role of geometrical non-uniformity in particular are far from clearly understood. Almost certainly, some of the features are a priori known, independent of the geometry (e.g. a more viscous displacing fluid provides a more efficient removal of the displaced fluid). Thus, an important question to address is when a geometrical non-uniformity becomes significant in governing the flow motion. In the light of the challenges and the question mentioned, some of the contributions of our model can be summarized as follows. (i) We take a first natural step to systematically study and understand the effects of a geometrical non-uniformity on displacement flows when various common features such as density ratios, viscosity ratios and shear-thinning and yield-stress parameters are present. (ii) Our semi-analytical model is useful in reducing the number of dimensionless parameters, making the study of these flows possible via fast computations. (iii) Our model helps to gain physical understanding of important laminar displacement flows by delivering regime classifications, interface behaviours, displacement front shapes, front velocities, displacement efficiency, etc. versus the dimensionless parameters. These results are valuable in providing guidance for future detailed experiments and computational fluid dynamic (CFD) computations on the flow regimes of interest. (iv) We focus mostly on the situations where a gentle convergence or divergence of the channel may have a significant effect on flow behaviours. This approach helps the identification of the mechanism by which we slowly deviate from the flow regimes previously observed and well studied for uniform channel flows. (v) Recent studies show the usefulness of small geometrical non-uniformities in controlling important displacement flow patterns, e.g. in HS flows (see Al-Housseiny et al. Reference Al-Housseiny, Tsai and Stone2012). In a related context, we would like to provide an understanding of the effects of geometrical non-uniformities on displacement flows for which the fluids coexist within the channel gap, while also adding the effects of buoyancy and non-Newtonian rheology. For the latter, considerable deviation from Newtonian flow behaviours may be expected. (vi) Simplified models such as ours are often developed to enable the stability analyses, permitting the prediction of flow regime transitions (see e.g. Alba, Taghavi & Frigaard (Reference Alba, Taghavi and Frigaard2013b ), who have recently extended the thin-film approach of Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009) into the weakly inertial range and performed a stability analysis).
The outline of the paper is as follows. In § 2, the geometry of the problem, the governing equations and the lubrication model are discussed. Section 3 presents the numerical procedure implemented. In § 4, interfacial behaviours are discussed and various flow regimes are classified for Newtonian displacement flows in a non-uniform channel. In § 5, the effects of shear-thinning and yield-stress properties in combination with the channel non-uniformity are discussed. Section 6 concludes the paper with a brief discussion and summary.
2 Flow geometry and governing equations
We consider in this problem a two-dimensional region between two plates, which are separated by distance
$\hat{D}_{0}+\tan {\it\alpha}\hat{x}\approx \hat{D}_{0}+{\it\alpha}\hat{x}$
(assuming
$|{\it\alpha}|\ll 1$
). The lower wall of the channel is oriented at an angle
${\it\beta}\approx {\rm\pi}/2$
with respect to the vertical. We take into account two mechanically stable types of displacement flows: a heavy fluid displacing a light fluid (1) in a slightly converging channel (for negative
${\it\alpha}$
) and (2) in a slightly diverging channel (for positive
${\it\alpha}$
). The downstream region between the two plates is initially filled with the lighter fluid (fluid
$L$
), which is displaced by a heavier fluid (fluid
$H$
). The latter is injected at a distance far away from the initial interface between the two fluids localized to
$\hat{x}=0$
. The mean imposed velocity is
$\hat{V}_{0}$
at
$\hat{x}=0$
. Cartesian coordinates
$(\hat{x},{\hat{y}})$
and the geometrical parameters are depicted in figure 1. We assume that the fluids are of generalized Newtonian type, for which the rheologies will be described later. The fluids studied are miscible but we consider the large-Péclet-number limit, the consequence of which is that no significant mixing occurs between the two fluids over the time scales of interest. We have made the Navier–Stokes equations dimensionless using the channel height at
$\hat{x}=0$
(i.e.
$\hat{D}_{0}$
) as length scale and
$\hat{V}_{0}$
as velocity scale. We have scaled time with
$\hat{D}_{0}/\hat{V}_{0}$
and pressure and stresses with
$\hat{{\it\mu}}_{H}\hat{V}_{0}/\hat{D}_{0}$
. We can therefore write the model equations as
$$\begin{eqnarray}\displaystyle & \displaystyle [1\pm At]Re[\boldsymbol{u}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}]=-\boldsymbol{{\rm\nabla}}p+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\tau}\pm \frac{Re}{Fr^{2}}\boldsymbol{e}_{g}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}=(u,v)$
denotes the velocity,
$p$
the pressure and
${\it\tau}$
the deviatoric stress. Also
$\boldsymbol{e}_{g}=(\cos {\it\beta},-\sin \,{\it\beta})$
is in directions
$(x,y)$
and the
$\pm$
refers to the heavy and light fluid layers, respectively. Note that we have subtracted the mean static pressure gradient from the pressure before scaling. The interface height is represented by
$y=h(x,t)$
. There are four dimensionless parameters that appear in (2.1). These are the inclination angle
${\it\beta}$
, the Atwood number defined as
$At=(\hat{{\it\rho}}_{H}-\hat{{\it\rho}}_{L})/(\hat{{\it\rho}}_{H}+\hat{{\it\rho}}_{L})$
, as well as the Reynolds number,
$Re$
, and the densimetric Froude number,
$Fr$
, which are defined as 
$$\begin{eqnarray}\displaystyle Re\equiv \frac{\hat{V}_{0}\hat{D}_{0}}{\hat{{\it\nu}}},\quad Fr\equiv \frac{\hat{V}_{0}}{\sqrt{At{\hat{g}}\hat{D}_{0}}}, & & \displaystyle\end{eqnarray}$$
where
$\hat{{\it\nu}}$
is defined using the mean density
$\hat{{\it\rho}}=(\hat{{\it\rho}}_{H}+\hat{{\it\rho}}_{L})/2$
and a viscosity scale derived from the rheological properties of pure fluid
$H$
. We are interested in flows with a small
$At$
, for which the flow is governed by the three dimensionless parameters
${\it\beta}$
,
$Re$
and
$Fr$
, a relevant dimensionless combination of which can be written as
$$\begin{eqnarray}\displaystyle {\it\chi}=\frac{2Re\cos {\it\beta}}{Fr^{2}}, & & \displaystyle\end{eqnarray}$$
which represents the balance of viscous stresses (due to the imposed flow) and axial buoyancy stresses (due to the density difference).
Figure 1. Schematic of the displacement flow geometry considered. Here
${\it\alpha}$
is negative so that the heavy fluid displaces the light one in a slightly converging channel. For a positive value of
${\it\alpha}$
, the channel would be slightly diverging. Note that
$|{\it\alpha}|\ll 1$
.
Before we proceed further, it should be clarified that we are not much interested in miscible displacement flows in microfluidic devices where transverse dimensions of flow geometry are extremely small, buoyant forces are negligible and dispersive regimes may be prevalent. We are instead motivated by laboratory-scale laminar displacement flow experiments and practical industrial-scale processes, which usually involve aqueous liquids with small density differences (e.g. 10 %) in geometries with a transverse dimension
$\hat{D}_{0}\sim 1$
cm and mean velocities
$\hat{V}_{0}\lesssim 10~\text{cm}~\text{s}^{-1}$
. These flows essentially belong to the category of high-
$Pe$
flows (typically
$Pe\geqslant 10^{6}$
) and they behave similarly to their immiscible analogues at infinite capillary number (i.e. vanishing surface tension). In the absence of instability, mixing and dispersion, the displacement flow interface remains relatively sharp over experimental time scales. For such flows, the Taylor-dispersion regime (Taylor Reference Taylor1953; Aris Reference Aris1956) is strictly found only for the dimensionless channel lengths
$\gg Pe$
. On the other hand, owing to the small density differences considered, the effect of
$At$
itself on the displacement flow is minor, although significant buoyancy effects (captured by
$Re/Fr^{2}$
) still exist. Consequently, our study of displacement flows of significant buoyancy effects at the non-dispersive high-
$Pe$
limit has practical relevance.
As usual, no-slip conditions at walls are satisfied. The channel is assumed to be sufficiently long in direction
$x$
. Thanks to our scaling, we can have in each cross-section of the channel
$$\begin{eqnarray}\displaystyle \int _{0}^{1+{\it\alpha}x}u\,\text{d}y=1. & & \displaystyle\end{eqnarray}$$
For buoyancy-dominated flows, the heavier fluid is expected to propagate at the bottom of the channel and the lighter fluid at the top. An interface, denoted by
$y=h(x,t)$
, separates the two pure fluids (see figure 1). We consider a two-layer flow so that the interface is single-valued. We assume in a usual fashion that velocity and stress are continuous across the interface, which is simply advected with the flow. Thus, the interface satisfies a kinematic condition.
2.1 Constitutive laws for Herschel–Bulkley fluids
We consider our fluids to be generalized Newtonian fluids, governed by a Herschel–Bulkley model, which comprises shear-thinning and yield-stress effects. This model also includes the Bingham model, as well as the power-law and Newtonian models. Constitutive laws for Herschel–Bulkley fluids are
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\gamma}}(\boldsymbol{u})=0\quad \Longleftrightarrow \quad {\it\tau}_{k}(\boldsymbol{u})\leqslant B_{k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{k,ij}(\boldsymbol{u})=\left[{\it\kappa}_{k}\dot{{\it\gamma}}^{n_{k}-1}(\boldsymbol{u})+\frac{B_{k}}{\dot{{\it\gamma}}(\boldsymbol{u})}\right]\dot{{\it\gamma}}_{ij}(\boldsymbol{u})\quad \Longleftrightarrow \quad {\it\tau}_{k}(\boldsymbol{u})>B_{k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{{\it\gamma}}_{ij}(\boldsymbol{u})=\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}, & & \displaystyle\end{eqnarray}$$
and the norms of these tensors,
$\dot{{\it\gamma}}(\boldsymbol{u})$
and
${\it\tau}_{k}(\boldsymbol{u})$
, are defined by
$$\begin{eqnarray}\displaystyle \dot{{\it\gamma}}(\boldsymbol{u})=\left[\frac{1}{2}\mathop{\sum }_{i,j=1}^{2}[\dot{{\it\gamma}}_{ij}(\boldsymbol{u})]^{2}\right]^{1/2},\quad {\it\tau}_{k}(\boldsymbol{u})=\left[\frac{1}{2}\mathop{\sum }_{i,j=1}^{2}[{\it\tau}_{k,ij}(\boldsymbol{u})]^{2}\right]^{1/2}. & & \displaystyle\end{eqnarray}$$
Here, subscripts
$k=H,L$
are used to distinguish the light and heavy fluids. Three dimensional parameters describe Herschel–Bulkley fluids, i.e. a fluid consistency
$\hat{{\it\kappa}}$
, a yield stress
$\hat{{\it\tau}}_{Y}$
and a power-law index
$n$
. Owing to our scaling, the parameter
${\it\kappa}_{H}=1$
and
${\it\kappa}_{L}$
is the viscosity ratio
$m$
defined as
$$\begin{eqnarray}\displaystyle m\equiv \frac{\hat{{\it\mu}}_{L}}{\hat{{\it\mu}}_{H}}=\frac{\hat{{\it\kappa}}_{L}[\hat{V}_{0}/\hat{D}_{0}]^{n_{L}-1}}{\hat{{\it\kappa}}_{H}[\hat{V}_{0}/\hat{D}_{0}]^{n_{H}-1}}, & & \displaystyle\end{eqnarray}$$
where
$\hat{{\it\mu}}_{L}$
is a viscosity scale for fluid
$L$
. Note that for two Newtonian fluids, we have
$\hat{{\it\mu}}_{k}=\hat{{\it\kappa}}_{k}$
. The Bingham numbers
$B_{k}$
are conventionally defined as
$$\begin{eqnarray}\displaystyle B_{k}\equiv \frac{\hat{{\it\tau}}_{k,Y}}{\hat{{\it\kappa}}_{H}[\hat{V}_{0}/\hat{D}_{0}]^{n_{H}}}. & & \displaystyle\end{eqnarray}$$
The scale for viscosity adopted in this paper is conventional and enables direct comparisons between our results and those of previous studies on displacement flows (e.g. Taghavi et al. Reference Taghavi, Seon, Martinez and Frigaard2009; Alba et al. Reference Alba, Taghavi and Frigaard2013b ). However, an alternative scaling approach is presented in appendix A.
2.2 Lubrication model
We are interested in a regime where the interface between the two fluids is elongated and where inertia is not dominant. The consideration of a near-horizontal angle has been due to the fact that we have previously found nearly viscous regimes at these angles in our experiment in uniform confined geometries (see e.g. Taghavi et al.
Reference Taghavi, Seon, Martinez and Frigaard2010, Reference Taghavi, Seon, Wielage-Burchard, Martinez and Frigaard2011, Reference Taghavi, Alba and Frigaard2012a
,Reference Taghavi, Alba, Seon, Wielage-Burchard, Martinez and Frigaard
c
). For example, the comparison between the lubrication model and experiments showed that, even at relatively significant Reynolds numbers, the displacement front velocities measured through experiments follow a viscous velocity scale and they can be approximated using the lubrication approach. Since in this paper we are dealing with a simple non-uniformity, i.e. a slightly converging or diverging channel, we may also expect here to find nearly viscous regimes even at relatively significant
$Re$
.
In order to scale our equations, we assume that the interface elongation is over a dimensionless length scale,
${\it\delta}^{-1}\gg 1$
, with a size of
$\hat{L}/\hat{D}_{0}$
, where
$\hat{L}$
is the characteristic spreading length. Aiming to derive a perturbation approximation, we rescale following standard methods using
${\it\delta}x=X$
,
${\it\delta}t=T$
,
${\it\delta}p=P$
and
$v={\it\delta}V$
to obtain the following reduced system of equations:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\it\delta}[1\pm At]Re\left[\frac{\partial u}{\partial T}+u\frac{\partial u}{\partial X}+V\frac{\partial u}{\partial y}\right]=-\frac{\partial P}{\partial X}+\frac{\partial }{\partial y}{\it\tau}_{Xy}\pm \frac{{\it\chi}}{2}+O({\it\delta}^{2}),\\ \displaystyle {\it\delta}^{3}[1\pm At]Re\left[\frac{\partial V}{\partial T}+u\frac{\partial V}{\partial X}+V\frac{\partial V}{\partial y}\right]=-\frac{\partial P}{\partial y}\mp {\it\delta}\frac{{\it\chi}}{2}\tan {\it\beta}+O({\it\delta}^{2}),\\ \displaystyle \frac{\partial u}{\partial X}+\frac{\partial V}{\partial y}=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
When
${\it\delta}\rightarrow 0$
with
$Re$
fixed, we can find
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\partial P}{\partial X}+\frac{\partial }{\partial y}{\it\tau}_{k,Xy}\pm \frac{{\it\chi}}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\partial P}{\partial y}\mp {\it\delta}\frac{{\it\chi}}{2}\tan {\it\beta}. & \displaystyle\end{eqnarray}$$
$(\hat{{\it\rho}}_{H}-\hat{{\it\rho}}_{L}){\hat{g}}\sin {\it\beta}\,\hat{D}_{0}$
, which act via the slope of the interface
${\sim}\hat{D}_{0}/\hat{L}$
. Considering that these stresses are balanced by viscous stresses gives 
$$\begin{eqnarray}\displaystyle {\it\delta}^{-1}=\frac{\hat{L}}{\hat{D}_{0}}=\frac{(\hat{{\it\rho}}_{H}-\hat{{\it\rho}}_{L}){\hat{g}}\sin {\it\beta}\,\hat{D}_{0}^{2}}{\hat{{\it\mu}}_{H}\hat{V}_{0}}={\it\chi}\tan {\it\beta} & & \displaystyle\end{eqnarray}$$
(see e.g. Alba et al.
Reference Alba, Taghavi and Frigaard2013b
). Therefore, the last term in (2.14) simplifies to
$1/2$
.
Assuming
${\it\delta}\rightarrow 0$
is of course a mathematical limit to apply the lubrication approximation at the leading order. In an experimental setting,
${\it\delta}$
needs to be small (meaning that
$\hat{L}/\hat{D}_{0}\gg 1$
) so that the model results can be compared with those of experiments. By rewriting the length scale versus the dimensionless parameters, we arrive at
${\it\delta}^{-1}=\hat{L}/\hat{D}_{0}=(Pe\sin {\it\beta}/ScFr^{2})\gg 1$
, where
$Sc=\hat{{\it\nu}}/\hat{D}_{m}$
is the Schmidt number, leading to
$Pe\sin {\it\beta}\gg ScFr^{2}$
. For a typical set of experimental parameters we have
$Pe\approx 10^{6}$
and
$ScFr^{2}\approx 2\times 10^{3}$
, satisfying the relation mentioned.
Another interpretation of the buoyancy parameter
${\it\chi}$
is a measure for the relative significance of the slope of the channel to that of the interface. For near-horizontal angles, the slopes of the interface and the channel are of tantamount importance. Therefore,
${\it\chi}$
may be regarded to be a parameter of order unity. For very large
${\it\chi}$
, the model may be still valid; however, a wrong scaling has been chosen. The case of very large
${\it\chi}$
can be interpreted as higher inclinations or smaller mean imposed velocities. In both cases,
$\hat{V}_{0}$
may not be the appropriate velocity scale. Instead, the appropriate velocity scale would be a characteristic viscous velocity that is obtained by balancing viscous and buoyant stresses. Therefore, the correct dimensionless parameter would be reduced to something like
$\cos {\it\beta}$
, as the channel slope has a dominant effect. In this work, we assume
${\it\chi}\geqslant 0$
, which may be interpreted as the slope of the channel being downhill in the direction of the flow.
Integrating (2.14) across both fluid layers delivers the pressure as
$$\begin{eqnarray}\displaystyle P(X,y,T)=\left\{\begin{array}{@{}ll@{}}\displaystyle P_{0}(X,T)+X\frac{{\it\chi}}{2}-\frac{y}{2}, & y\in [0,h],\\ \displaystyle P_{0}(X,T)+X\frac{{\it\chi}}{2}+\frac{(y-2h)}{2}, & y\in [h,1+{\it\alpha}X],\end{array}\right. & & \displaystyle\end{eqnarray}$$
where
$P_{0}(X,T)$
is defined by
$P_{0}(X,T)=P(X,0,T)-X({\it\chi}/2)$
. On substituting into (2.13), we arrive at
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\partial P_{0}}{\partial X}+\frac{\partial }{\partial y}{\it\tau}_{H,Xy},\quad y\in (0,h), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\partial P_{0}}{\partial X}+\frac{\partial }{\partial y}{\it\tau}_{L,Xy}-{\it\chi}+\frac{\partial h}{\partial X},\quad y\in (h,1+{\it\alpha}X). & \displaystyle\end{eqnarray}$$
The leading-order strain-rate component in the lubrication approximation is
$\dot{{\it\gamma}}_{Xy}=\partial u/\partial y$
. The leading-order shear stress,
${\it\tau}_{k,Xy}$
, is also defined in terms of
$\dot{{\it\gamma}}_{Xy}$
using the following leading-order constitutive laws for pure fluid
$k$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u}{\partial y}=0\quad \Longleftrightarrow \quad |{\it\tau}_{k,Xy}|\leqslant B_{k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{k,Xy}=\left[{\it\kappa}_{k}\left|\frac{\partial u}{\partial y}\right|^{n_{k}-1}+\frac{B_{k}}{\left|\displaystyle \frac{\partial u}{\partial y}\right|}\right]\frac{\partial u}{\partial y}\quad \Longleftrightarrow \quad |{\it\tau}_{k,Xy}|>B_{k}. & \displaystyle\end{eqnarray}$$
Equations (2.17) and (2.18) define an elliptic problem for
$u(y)$
, for given
$h$
and
$\partial h/\partial X$
. For
$u(y)$
, the boundary conditions are
$u=0$
at the lower wall (
$y=0$
) and at the upper wall (
$y=1+{\it\alpha}X$
). At the interface between the two fluids (
$y=h$
) the velocity is continuous. The stress continuity can also be represented by
${\it\tau}_{H,Xy}={\it\tau}_{L,Xy}$
. The conditions stated allow for the determination of
$u$
for a given
$\partial P_{0}/\partial X$
. Equation (2.5) has to be satisfied as the additional constraint to determine the pressure gradient. The interface is advected via a kinematic condition furnishing the dependence of the velocity on space and time:
$$\begin{eqnarray}\displaystyle \frac{\partial h}{\partial T}+u\frac{\partial h}{\partial X}=V. & & \displaystyle\end{eqnarray}$$
Combining the incompressibility condition with the kinematic condition results in
$$\begin{eqnarray}\displaystyle \frac{\partial h}{\partial T}+\frac{\partial q}{\partial X}=0, & & \displaystyle\end{eqnarray}$$
where
$q$
is defined through
$$\begin{eqnarray}\displaystyle q=\int _{0}^{h}u\,\text{d}y. & & \displaystyle\end{eqnarray}$$
In the rest of our paper we concentrate on providing the solutions to the system (2.22) and (2.23).
Regarding boundary conditions, we assume that the channel is full of pure fluid
$L$
and fluid
$H$
at its two ends. Concerning initial conditions, we must consider an initial profile compatible with the far-field conditions, leading to
$$\begin{eqnarray}\displaystyle h(X,0)\rightarrow 1+{\it\alpha}X-H(X)(1+{\it\alpha}X). & & \displaystyle\end{eqnarray}$$
Here
$H(X)$
is the usual Heaviside function, used to ensure that the initial variation in
$h$
localized to
$X=0$
is sharp. Note that we always limit the length of the channel so that the channel upper wall never contacts the lower one. This is intuitively required to conform to the inlet and exit conditions. For a diverging channel (positive
${\it\alpha}$
) we limit the length to
$X\in (-{\it\alpha}^{-1},\infty )$
and for a converging channel (negative
${\it\alpha}$
) to
$X\in (-\infty ,-{\it\alpha}^{-1})$
. In the following sections, for simplicity we may refer to
$X$
as distance, providing a measure of a distance from
$X=0$
, the spatial location of the initial sharp variation in
$h$
.
2.3 Computing the flux function
We assume that, for physically sensible dimensionless parameters, there is always a velocity solution, from which the flux can be computed. In this section we explain how this can be done in practice. First, note that, for fixed
$h$
and
$\partial h/\partial X$
, equations (2.17) and (2.18) may be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial y}{\it\tau}_{H,Xy}=\frac{\partial P_{0}}{\partial X},\quad y\in (0,h), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial y}{\it\tau}_{L,Xy}={\it\chi}-h_{X}+\frac{\partial P_{0}}{\partial X},\quad y\in (h,1+{\it\alpha}X). & \displaystyle\end{eqnarray}$$
$y$
, in each pure fluid layer. For Newtonian displacement flow, the analytical solution may be easily found as 
$$\begin{eqnarray}\displaystyle q(h,{\it\chi},h_{X},m,{\it\alpha}X)=\frac{\displaystyle \mathop{\sum }_{i=0}^{4}f_{i}({\it\alpha}X)^{i}}{\displaystyle \mathop{\sum }_{i=0}^{4}g_{i}({\it\alpha}X)^{i}}, & & \displaystyle\end{eqnarray}$$
where the functions
$f_{i}$
are
$$\begin{eqnarray}\displaystyle f_{0} & = & \displaystyle {\textstyle \frac{1}{3}}h^{2} (9m-4{\it\chi}h^{4}+h^{2}m{\it\chi}-3h^{3}m{\it\chi}+3h^{4}m{\it\chi}-h^{2}mh_{X}+3h^{3}mh_{X}-3h^{4}mh_{X}\nonumber\\ \displaystyle & & \displaystyle +\,h^{5}mh_{X}-hh_{X}+4h^{2}h_{X}+4h_{X}h^{4}-3h^{2}m+h^{5}{\it\chi}+h{\it\chi}-4h^{2}{\it\chi}-h^{5}m{\it\chi}+6{\it\chi}h^{3}\nonumber\\ \displaystyle & & \displaystyle -\,6h_{X}h^{3}+3h^{2}m^{2}-6hm-h^{5}h_{X} ),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle f_{1} & = & \displaystyle {\textstyle \frac{1}{3}}h^{2} (12{\it\chi}h^{3}+4h_{X}h^{4}+3h^{2}m{\it\chi}-6h^{3}m{\it\chi}+12h^{2}h_{X}-3h^{4}mh_{X}-6hm-3h^{2}mh_{X}\nonumber\\ \displaystyle & & \displaystyle +\,3h^{4}m{\it\chi}-12h^{2}{\it\chi}+4h{\it\chi}-12h_{X}h^{3}-4hh_{X}+6h^{3}mh_{X}+18m-4{\it\chi}h^{4} ),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle f_{2} & = & \displaystyle {\textstyle \frac{1}{3}}h^{2} (3h^{3}mh_{X}-12h^{2}{\it\chi}+12h^{2}h_{X}+9m+6h{\it\chi}-3h^{3}m{\it\chi}-3h^{2}mh_{X}+6{\it\chi}h^{3}\nonumber\\ \displaystyle & & \displaystyle +\,3h^{2}m{\it\chi}-6hh_{X}-6h_{X}h^{3} ),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle f_{3} & = & \displaystyle {\textstyle \frac{1}{3}}h^{2}(-h^{2}mh_{X}+h^{2}m{\it\chi}+4h^{2}h_{X}-4h^{2}{\it\chi}-4hh_{X}+4h{\it\chi}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle f_{4} & = & \displaystyle {\textstyle \frac{1}{3}}h^{2}(h{\it\chi}-hh_{X}).\end{eqnarray}$$
$g_{i}$
are obtained as 
$$\begin{eqnarray}\displaystyle g_{0} & = & \displaystyle 1-4h+h^{4}+h^{4}m^{2}-2h^{4}m+4hm-4h^{3}-6h^{2}m+4h^{3}m+6h^{2},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle g_{1} & = & \displaystyle 4-12h+4h^{3}m+12hm-12h^{2}m+12h^{2}-4h^{3},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle g_{2} & = & \displaystyle 12hm+6-12h+6h^{2}-6h^{2}m,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle g_{3} & = & \displaystyle 4hm-4h+4,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle g_{4} & = & \displaystyle 1.\end{eqnarray}$$
${\it\alpha}=0$
, we find
$q=f_{0}/g_{0}$
, which is the same as the flux function for a uniform channel flow, for which we have successfully compared our result with that given in Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009).
Figure 2 shows examples of contours of
$q$
versus
$h$
and
$X$
for various values of
${\it\alpha}$
and
$m$
, for Newtonian displacements. Comparing figures 2(a–c), for a fixed interface height in a uniform channel, we observe that the flux increases for larger
$m$
. On the other hand, for a fixed flux, the interface height is larger for smaller
$m$
. This implies that, even before performing numerical simulations, we may expect to find larger interface heights for smaller
$m$
, meaning that displacing fluid is relatively more successful in pushing the displaced fluid out of the channel geometry. The effect of the viscosity ratio for a converging channel flow is similar (see figure 2
d–f), although at larger
$X$
the interface height varies over such a narrow channel that the flux value becomes very sensitive to the exact value of the interface height. For flows occurring in a diverging channel (see figure 2
g–i), the flux can become larger than
$1$
for particular values of interface heights at certain distances (e.g. at
$X\gtrapprox 20$
in figure 2
g). This implies that, for these interface heights, there has to exist a local flow reversal of the displaced fluid against the direction of imposed flow to compensate for the large flux of the displacing fluid.
Figure 2. Examples of contours of
$q$
versus
$h$
and
$X$
for
${\it\chi}=20$
and
$h_{X}=0$
, for Newtonian displacements. In each row,
${\it\alpha}$
is fixed: 0 (a–c),
$-0.01$
(d–f), 0.01 (g–i). In each column,
$m$
is fixed: 0.1 (a,d,g), 1 (b,e,h), 10 (c,f,i).
The analytical form of the Newtonian flux function may be used as a test solution for the more general problem, to which we now turn. The wall shear stresses in fluids
$H$
and
$L$
are denoted by
${\it\tau}_{H}$
and
${\it\tau}_{L}$
, respectively. These can be defined versus the pressure gradient
$\partial P_{0}/\partial X$
and interfacial stress
${\it\tau}_{i}$
following
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{H}={\it\tau}_{i}-h\frac{\partial P_{0}}{\partial X}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{L}={\it\tau}_{i}+(1+{\it\alpha}X-h)\left({\it\chi}-h_{X}+\frac{\partial P_{0}}{\partial X}\right). & \displaystyle\end{eqnarray}$$
${\it\tau}_{i}$
,
${\it\tau}_{H}$
and
${\it\tau}_{L}$
in each layer are 
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{H,Xy}(y)={\it\tau}_{H}\left(1-\frac{y}{h}\right)+{\it\tau}_{i}\frac{y}{h}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{L,Xy}(y)={\it\tau}_{L}\frac{h-y}{h-1-{\it\alpha}X}+{\it\tau}_{i}\frac{1+{\it\alpha}X-y}{1+{\it\alpha}X-h}. & \displaystyle\end{eqnarray}$$
$\partial u/\partial y$
can be defined at each point using the constitutive laws governing the displacing and displaced fluids. It then becomes possible to integrate
$\partial u/\partial y$
away from the walls at
$y=0$
and
$y=1+{\it\alpha}X$
, where the no-slip conditions are applied, towards the interface. For given
${\it\tau}_{H}$
,
${\it\tau}_{L}$
and
${\it\tau}_{i}$
, and the rheological parameters of each fluid, two interfacial velocities can be found as 
$$\begin{eqnarray}\displaystyle & \displaystyle u_{i}(h^{-})=\int _{0}^{h}\frac{u(y;{\it\tau}_{H},{\it\tau}_{i})}{\partial y}\,\text{d}y, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{i}(h^{+})=\int _{1+{\it\alpha}X}^{h}\frac{u(y;{\it\tau}_{L},{\it\tau}_{i})}{\partial y}\,\text{d}y. & \displaystyle\end{eqnarray}$$
$({\it\tau}_{H},{\it\tau}_{L})$
. However, it is possible to simply iterate on
${\it\tau}_{i}$
to reach 
$$\begin{eqnarray}\displaystyle {\rm\Delta}u_{i}({\it\tau}_{i})\equiv u_{i}(h^{-})-u_{i}(h^{+})=0. & & \displaystyle\end{eqnarray}$$
For the case of a uniform channel (i.e.
${\it\alpha}=0$
), the detailed procedure is given in appendix A of Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009) and we follow the same approach, albeit for boundary conditions of a non-uniform channel. After some algebra, it is indeed possible to provide lengthy analytical expressions for the two interfacial velocities (versus
${\it\tau}_{i}$
,
${\it\tau}_{H}$
and
${\it\tau}_{L}$
) and, subsequently, for the flux functions in each fluid layer, for a non-uniform channel flow. For brevity, we do not present these here.
Examples of contours of
$q$
versus
$h$
and
$X$
for non-Newtonian displacements are shown in figure 3. It can be seen that the non-Newtonian rheology can significantly affect the variation of the flux in both converging and diverging channels. Regarding a converging channel flow (see figure 3
a,b), the shape of the contours is modified depending on which fluid (displacing or displaced) is shear-thinning. In the diverging situation (see figure 3
e,f), the maximum flux is seen at large distances at an intermediate interface height, the value of which depends on which fluid is shear-thinning. Yield-stress effects can modify the flux contours. A particularly interesting example is when a Newtonian fluid displaces a yield-stress fluid in a diverging channel (i.e. figure 3
h), where one observes a large contour region with the flux almost equal to 1. The height at which the flux reaches unity seems to be only slowly increasing over a large length of the diverging channel.
Figure 3. Examples of contours of
$q$
versus
$h$
and
$X$
for
$m=1$
,
${\it\chi}=20$
and
$h_{X}=0$
, for non-Newtonian displacements. In each row,
${\it\alpha}$
is fixed:
$-0.01$
(a–d), 0.01 (e–h). In each column,
$[n_{H},n_{L},B_{H},B_{L}]$
$=$
$[0.25,1,0,0]$
(a,e),
$[1,0.25,0,0]$
(b,f),
$[1,1,20,0]$
(c,g),
$[1,1,0,20]$
(d,h).
2.4 Remarks
Before proceeding, we comment on a few points for clarification and justification of the direction taken for the rest of the paper.
-
(i) For a uniform channel boundary condition, Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009) showed that there are no steady travelling wave solutions to (2.22). This fact has been implicitly discussed in Seon et al. (Reference Seon, Znaien, Salin, Hulin, Hinch and Perrin2007a ) for an exchange flow of Newtonian fluids with identical viscosities. These arguments can be rigorously extended to a slightly non-uniform channel flow. However, we proceed, without a proof, to consider that (2.22) has no steady travelling wave solutions with non-uniform channel boundary conditions. This implies that, under the lubrication assumption, for no combination of the parameters involved, a complete removal of the displaced fluid by the displacing fluid can be achieved. As a consequence, in general, numerical simulations of the interface evolution in time must be conducted to survey as many parameters as possible, in order to identify which displacement is ‘better’ or ‘worse’. The procedure for numerical simulations is given in § 3.
-
(ii) One may be interested by the main features of displacement flows that can be studied using a lubrication model, e.g. long- and short-time behaviours of the interface (cf. Taghavi et al. Reference Taghavi, Seon, Martinez and Frigaard2009), stationary interface flows (cf. Taghavi et al. Reference Taghavi, Seon, Wielage-Burchard, Martinez and Frigaard2011), front height and speeds (cf. Seon et al. Reference Seon, Znaien, Salin, Hulin, Hinch and Perrin2007a ; Taghavi et al. Reference Taghavi, Alba, Seon, Wielage-Burchard, Martinez and Frigaard2012c ) and displacement efficiencies. One approach to study these features is directly solving the kinematic equation for a given set of parameters to simulate the spatial and temporal evolution of the interface. Nevertheless, as discussed in § 4.2, some of the chief features of the interface may be approximated using the long-time limit assumption and without fully solving the equations.
-
(iii) There are seven dimensionless parameters that govern the problem that we study, i.e.
${\it\chi}$
,
$m$
,
$B_{k}$
,
$n_{k}$
and
${\it\alpha}X$
. This relatively large number of parameters makes it hard to deliver quantitative predictions for all the possible flow features (e.g. those mentioned above). Therefore, in order to provide essential understanding of displacement flows in non-uniform channels, we will consider the Newtonian and non-Newtonian fluids separately, focusing on varying
$m,{\it\chi}$
and the rheological parameters, respectively. -
(iv) This paper is focused on density-unstable displacements (i.e. a heavy fluid displacing a light fluid). Although, with almost no effort, the model presented in this paper can also be exploited to obtain the results for density-stable displacements (i.e. a light fluid displacing a heavy fluid), this has not been done, for the following reasons. First, instead of adding results, we would like to keep the paper focused on the effects of the flow geometry. Second, while there are good comparisons between density-unstable displacement experiments and the associated lubrication model (see e.g. Taghavi et al. Reference Taghavi, Alba, Seon, Wielage-Burchard, Martinez and Frigaard2012c ), comparisons between the lubrication model and density-stable experiments are limited: unlike the model, experiments show that the trailing front is not pinned to the wall and it moves downstream (see e.g. Alba et al. Reference Alba, Taghavi and Frigaard2012).
3 Numerical procedure
For numerically solving our equations, we suppose fully developed flows of pure heavy and light fluids at the two ends of the channel. We have observed that the interface propagation at
$T\sim O(1)$
is insensitive to reasonably sharp initial conditions, taken as a linear function of
$X$
: typically
$h(X,T=0)=-X+0.5$
. The kinematic condition (2.22) is discretized in conservative form, second-order in space and first-order explicitly in time. It is integrated using a Van Leer flux limiter scheme (see e.g. Yee, Warming & Harten Reference Yee, Warming and Harten1985) for shock capturing.
Figure 4. Dependence of the interface evolution on mesh size in a converging channel: (a) a Newtonian fluid displacing another Newtonian fluid; (b) a Newtonian fluid displacing a yield-stress fluid with
$B_{L}=5$
and
$n_{L}=1$
; (c) a Newtonian fluid displacing a shear-thinning fluid with
$n_{L}=0.5$
. The other parameters used are
${\it\alpha}=-0.01$
,
$m=1$
,
${\it\chi}=-20$
at
$T=35$
, with mesh steps
$\text{d}X=0.02$
(dash-dotted line),
$\text{d}X=0.05$
(dashed line) and
$\text{d}X=0.1$
(solid line). The time step is fixed to
$dT=0.001$
in all panels. The insets are zoomed at the fronts.
To explore the performance of the numerical procedure, figure 4 examines convergence with varying spatial mesh step
$\text{d}X=0.02,0.5,0.1$
. The inclined thick lines in figure 4 (and in the rest of the figures in this paper) indicate the position of the inclined upper wall. The interface evolutions observed result in an advancing displacement front, stretching out the upper part of the interface. We have selected one Newtonian and two non-Newtonian displacements. Figure 4 illustrates that the interface propagation does not depend on the choice of the mesh step as long as
$\text{d}X$
is selected to be sufficiently small. However, the mesh step has little effects on the front of the interface; see the insets in figure 4. We have benchmarked our results against those for uniform channel flows of Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009, Reference Taghavi, Seon, Martinez and Frigaard2010). Around 300 simulations have been performed, from which the main findings are presented in this paper.
4 Newtonian fluids
We start by exploring Newtonian fluid displacements, mainly for three reasons. First, some of the major qualitative behaviours are observed in Newtonian fluid displacements. Second, having fewer dimensionless parameters greatly simplifies the analysis of these fluids. Finally, numerical solution is markedly faster due to the analytical expression obtained for the flux function.
In this work, for simplicity, we fix
${\it\alpha}$
to two small values, i.e.
${\it\alpha}=-0.01$
(converging) and
${\it\alpha}=0.01$
(diverging), and compare the results against those of uniform channel flows (
${\it\alpha}=0$
).
4.1 Examples of typical qualitative behaviour
Examples of displacements of Newtonian fluids for a fixed
${\it\chi}$
and various values of
$m$
are plotted in figure 5. The results are shown at
${\it\chi}=0$
, which is a representative case for either of the following situations: (i) the channel is strictly horizontal (
${\it\beta}={\rm\pi}/2$
); (ii) the channel is slightly inclined but
${\it\chi}$
is very small, meaning that viscous forces are relatively very large. Our simulations (not shown here) confirm that the results for
${\it\chi}=0$
,
$0.1$
and
$0.5$
are in fact identical. The length of domain is the same in all of the panels in figure 5, so that a direct visual comparison between the results in different geometries becomes possible. In uniform channel displacement flows in figure 5(a–c), the initially steep interface elongates as time grows. Our simulations demonstrate that the interface between the two fluids develops into a slumping profile, for which the solution generally has two segments. The first segment is an advancing leading front towards the bottom of the channel that has a constant shape and that moves at constant speed. At the top of the interface there is a second region, where the interface is stretched. This has a trailing front towards the top of the channel, which for the cases shown in figure 5 appears to be pinned to the upper wall. For larger
$m$
, the height of the interface is relatively smaller and its front speed is larger; therefore, the interface advances relatively farther over the same period of time. When the channel is converging, figure 5(d–f) show similarities but also differences compared with the uniform channel flow results. First of all, the interface analogously has the two segments as explained before. Nonetheless, as the interface evolves, its height decreases while its front speed increases. The effect of
$m$
is expectedly the same as observed for the uniform channel flow. Figure 5(g–i) suggest that the interfacial behaviour for the diverging channel displacement flow is the opposite of that for the converging one. The front height increases whereas the front speed decreases as the interface evolves. The dashed lines in each panel of figure 5 depict the front height predictions at long times using (4.4) (discussed in § 4.2), showing very good agreement. Having computed front heights from (4.4), it becomes possible to predict front speeds matching the simulation results as expected (although the results are not shown).
Figure 5. Examples of Newtonian displacements for
${\it\chi}=0$
and
$T=0,3,\ldots ,30$
. In each row,
${\it\alpha}$
is fixed: 0 (a–c),
$-0.01$
(d–f), 0.01 (g–i). In each column,
$m$
is fixed: 0.1 (a,d,g), 1 (b,e,h), 10 (c,f,i). The dashed lines illustrate the front height predictions of (4.4) for long times.
Figure 6. Interface evolution in time,
$T=0,3,\ldots ,27,30$
, for
$m=1$
. In each row,
${\it\alpha}$
is fixed: 0 (a–c),
$-$
0.01 (d–f), 0.01 (g–i). In each column,
${\it\chi}$
is fixed: 20 (a,d,g), 70 (b,e,h), 100 (c,f,i). The dashed lines illustrate the front height predictions of (4.4) for long times. For panel (f) only, a very small artificial diffusion has been added to stabilize the numerical method.
Figure 5 presented the results for fixed
${\it\chi}=0$
, which exemplifies cases of very large viscous stresses compared with buoyant stresses in the direction of the flow, e.g. large imposed flows and vanishingly small density differences. Figure 6 shows, on the other hand, the effect of increasing
${\it\chi}$
(buoyancy) for displacement flows in uniform and non-uniform channels. For relatively small
${\it\chi}$
, in figure 6(a,d,g) (for
${\it\chi}=20$
), the existence of buoyancy modifies the interface shape. However, compared to figure 5, the displacement patterns generally seem similar. As can be seen, for large
${\it\chi}$
in figure 6(c,f,i), the buoyancy force is so strong that it creates a motion of the displaced fluid against the flow direction. The trailing front is no longer pinned to the upper wall and it starts to move backwards. It is interesting to note that the onset of the appearance of the back-flow depends on the buoyancy parameter (
${\it\chi}$
) but not on the geometry parameter (
${\it\alpha}$
). For example, figure 6(b,e,h) display nearly stationary interfaces for all values of
${\it\alpha}$
at
${\it\chi}=70$
, which is approximately the critical
${\it\chi}$
for the onset of the back-flow. In addition, figure 6(c,f,i) illustrate that there exists a back-flow for all values of
${\it\alpha}$
at
${\it\chi}=100$
. Despite these similarities, in § 4.4 we will discuss a fundamental difference between displacements at large
${\it\chi}$
in converging and diverging channels by demonstrating that, for any given parameter set, there is a distance at which the back-flow stops only for the diverging channel flow (i.e. positive
${\it\alpha}$
).
The (red) dashed lines superimposed in figure 6 depict the front height predictions from (4.4) for long times. The predictions for the uniform channel flows in figure 6(a–c) are nearly perfect. The agreement of the predicted front heights and the simulations of converging channel flows are reasonable, although a deviation is noticeable as time grows. However, the results for the diverging channel flows in figure 6(g–i) tell a different story. While at
${\it\chi}=20$
, the agreement is not bad, at
${\it\chi}=70$
and
${\it\chi}=100$
, the predictions and simulation results do not agree and follow opposite paths, e.g. the predicted front heights increase with
$X$
whereas the simulation front heights decrease. Therefore, owing to the crude and approximate nature of the approach given in § 4.2, it cannot provide accurate predictions of the front height for diverging channel flows at large
${\it\chi}$
. The decrease in the front heights observed in figure 6(h,i) is actually a contour-intuitive feature, which we now attempt to explain. For large
${\it\chi}$
, the buoyancy is stronger. As time grows, the leading front penetrates into the displaced fluid so that it advances through a wider channel. Therefore, it locally experiences progressively larger buoyancy forces (proportional to the thickness of the channel) at larger
$X$
, due to which the front needs to move faster. The flux constraint instead restricts the fast front to have a shorter height. In fact, it is noteworthy that, among all the results presented in figure 6, the diverging channel flow front is the one that travels fastest and farthest!
We conclude this section with a few comments. Firstly, regarding the interface, one may anticipate to observe diffusive behaviours (concerning the interface slope) at short times. However, our investigation both analytically and numerically shows that the flow behaviour in a non-uniform channel at short times recovers that in a uniform channel. Therefore, the rest of our study will be concentrated on the longer-time behaviour of the interface. Secondly, long times in terms of our lubrication model are when
$T>1$
but we also typically have
$T<100$
, considering the typical length of the channel. This long time is much smaller than the time scale over which molecular diffusion dominates the flow. It can be shown that smearing of the interface and effective diffusion across the channel are observed at
$t_{diffusion}\sim Pe$
, eventually approaching the dispersive limit (see also Chen & Meiburg Reference Chen and Meiburg1996; Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Rakotomalala et al.
Reference Rakotomalala, Salin and Watzky1997; Yang & Yortsos Reference Yang and Yortsos1997). This time can be rewritten as
$T_{diffusion}\sim Fr^{2}Sc$
using our lubrication rescaling and assuming
$\sin {\it\beta}\approx 1$
. For typical physical parameters, we find that
$T_{diffusion}\sim O(10^{3}{-}10^{5})$
, which is much longer than the long times that we consider in the lubrication model results. Thirdly, the brief overview of the typical, qualitative examples of displacement flows showed that two main flow regimes may be distinguished thus far for both uniform and non-uniform channel flows through our lubrication model: a no-back-flow regime (for small values of
${\it\chi}$
), discussed further in § 4.3, and a sustained back-flow regime (for relatively large values of
${\it\chi}$
), discussed in § 4.4. The buoyancy force is weak in the no-back-flow regime while it is fairly strong in the sustained back-flow regime, causing the displaced fluid to move backwards in the opposite direction of the mean imposed flow. The transition between the two regimes, quantified by a critical value of
${\it\chi}$
for a given set of parameters, is marked by a stationary interface flow regime discussed further in § 4.4. Finally, the results for no-back-flow and sustained back-flow regimes are fundamentally different. For example, for the no-back-flow regime, § 4.3 reveals that it is possible to find a similarity form to superimpose the interface evolutions at very long times for both converging and diverging channel geometries, i.e. for both negative and positive
${\it\alpha}$
. However, such a superposition of results cannot be attained for the sustained back-flow regime in a non-uniform channel.
4.2 Behaviour of the interface at long times
Since we have directed our focus on the behaviour of the interface at long times, it is natural to attempt to directly compute this long-time behaviour (using some assumptions), in place of computations. We have already observed that a typical interface has two segments. Looking at the upper stretched region, we have seen that at long times the slope of the interface,
$h_{X}$
, is small, which we assume to be zero as
$T\gg 1$
. Thus, the interface evolution at long times in the stretched region should be governed approximately by the following hyperbolic (rather than parabolic) equation:
$$\begin{eqnarray}\displaystyle \frac{\partial h}{\partial T}+\frac{\partial \tilde{q}(h,{\it\alpha}X)}{\partial X}=0, & & \displaystyle\end{eqnarray}$$
where
$\tilde{q}=q_{h_{X}=0}$
. The interface is elongated between the location of the leading front,
$X_{f}$
, and that of the trailing front,
$X_{t}$
. The interface height between these two values varies in the interval of
$h\in [0,1+{\it\alpha}X_{t}]$
. At long times, as the total area behind the interface must conserve mass, we can write
$$\begin{eqnarray}\displaystyle \int _{0}^{1+{\it\alpha}X_{t}}\frac{\partial \tilde{q}(h,{\it\alpha}X)}{\partial h}\,\text{d}h=1. & & \displaystyle\end{eqnarray}$$
Additionally, introducing
$V_{f}$
and
$h_{f}$
, respectively, as the interface speed and the leading front height at
$X_{f}$
, we can rewrite
$$\begin{eqnarray}\displaystyle \int _{0}^{1+{\it\alpha}X_{t}}\frac{\partial \tilde{q}(h,{\it\alpha}X)}{\partial h}\,\text{d}h=V_{f}h_{f}+\int _{h_{f}}^{1+{\it\alpha}X_{t}}\frac{\partial \tilde{q}(h,{\it\alpha}X)}{\partial h}\,\text{d}h=1. & & \displaystyle\end{eqnarray}$$
Using a crude estimation for the leading front velocity as
$V_{f}\approx (\partial \tilde{q}/\partial h)(h_{f},{\it\alpha}X_{f})$
, we arrive at
$$\begin{eqnarray}\displaystyle \tilde{q}(h_{f},{\it\alpha}X_{f})\approx h_{f}\frac{\partial \tilde{q}}{\partial h}(h_{f},{\it\alpha}X_{f}). & & \displaystyle\end{eqnarray}$$
This approximate equation can now be used to estimate the front height
$h_{f}$
, for which the solution can be directly obtained through the use of the conventional equal areas rule. When the approximate value of the front height is obtained, the front velocity can be calculated straightforwardly. As discussed earlier, the dashed lines in figures 5 and 6 show the results of the front height predictions at long times through solving (4.4).
4.3 No-back-flow regime
The no-back-flow regime occurs at relatively smaller values of
${\it\chi}$
, i.e. for smaller ratios of buoyant to viscous stresses. For small
${\it\chi}$
, we have observed that the interface evolves quickly into a shape consisting of a front region and a stretched region. For
${\it\alpha}=0$
, the front region has a constant height and advances at steady speed in the lower part of the channel. For
${\it\alpha}<0$
(
${\it\alpha}>0$
) the front region shortens (expands) and accelerates (decelerates) as the flow develops. Regarding the stretched region for all values of
${\it\alpha}$
, the interface is continuously extended as time increases. Perhaps interesting is that there exists a similarity form for the interface evolution for both uniform and non-uniform channel flows, where the interfaces at long times can be superimposed. The longer-time profiles of
$h/(1+{\it\alpha}X)$
can be conveniently plotted against
$(X+{\it\alpha}X^{2}/2)/T$
, in which the interface profiles collapse onto a single similarity profile as
$T\gg 1$
(see figure 7). It is easy to obtain this similarity form. Postulating that the varying front speed for a non-uniform channel flow (
$V_{f}$
) may be related to the steady front speed for a corresponding uniform channel flow (
$V_{f}^{\ast }$
) through the relation
$(1+{\it\alpha}X_{f})V_{f}\approx V_{f}^{\ast }$
, we can replace
$V_{f}$
with
$\text{d}X_{f}/\text{d}T$
and integrate to find the similarity variable as
$(X+{\it\alpha}X^{2}/2)/T$
, where
$X_{f}$
has been substituted by
$X$
. On the other hand, we postulate that the front height at long times in a non-uniform channel flow may be found through the relation
$h_{f}/(1+{\it\alpha}X_{f})\approx h_{f}^{\ast }$
, where
$h_{f}^{\ast }$
is the constant interface height in a corresponding uniform channel flow. Thus, we can generalize and obtain that the interface height must be scaled as
$h/(1+{\it\alpha}X)$
to furnish the superposition.
At relatively long times, the following relation between the leading front speed and the location of the leading front can be straightforwardly obtained:
$$\begin{eqnarray}\displaystyle V_{f}\approx \frac{1}{T}\left(\frac{X_{f}}{(1+{\it\alpha}X_{f})}+\frac{{\it\alpha}{X_{f}}^{2}}{2(1+{\it\alpha}X_{f})}\right). & & \displaystyle\end{eqnarray}$$
For a converging channel flow, this relation shows, at small
$X_{f}$
, that the leading front speed may be approximated by
$X_{f}/T$
, which is the same as the constant leading front speed in a uniform channel. This may lead to the interpretation that the displacement front does not feel the existence of the wall non-uniformity at short distances. However, as the channel narrows down and the front advances, the front speed increases significantly and essentially approaches infinity as
$X_{f}\rightarrow -1/{\it\alpha}$
. We define a critical transition distance at which the front speed becomes twice as large as the speed at short distances as a breakthrough distance where the effect of the non-uniformity starts to become significantly felt by the leading front. This transition distance can be found as
$X_{trans}^{{\it\alpha}<0}\sim -2/3{\it\alpha}$
. On the other hand, for a diverging channel flow, the front speed varies between approximately
$X_{f}/T$
at short distances and
$X_{f}/2T$
at very long distances (as
$X_{f}\rightarrow \infty$
). In this case, we may define a breakthrough transition distance as
$X_{trans}^{{\it\alpha}>0}\sim 1/{\it\alpha}$
at which the front speed drops below the average of the speeds at short and long distances. In fact, our simulations (not shown here) confirm that the behaviour of the flow at the breakthrough distances starts to change. Nevertheless, we acknowledge that no single measure or definition might be universal for breakthrough distances of such displacement flows.
Figure 7. Superposition of the interface heights in Newtonian displacements for
${\it\chi}=0$
and
$T=3,\ldots ,30$
by plotting
$h/(1+{\it\alpha}X)$
as a function of
$(X+{\it\alpha}X^{2}/2)/T$
for
${\it\alpha}=-0.01$
(solid line),
${\it\alpha}=0$
(dashed line) and
${\it\alpha}=0.01$
(dash-dotted line):
$m=0.1$
(a), 1 (b) and 10 (c). Note that the same parameters as in figure 5 have been used. The insets are for
$T=0.3,0.6,\ldots ,3$
.
Let us take another look at figure 7 showing a nearly perfect superposition of the interface heights versus the similarity variable for the same parameters as in figure 5 for all the values of
${\it\alpha}$
studied. This finding has an important implication: for the no-back-flow regime it is possible to extract the behaviour of the interface in a converging or diverging channel flow through the solution readily obtained for the corresponding uniform channel flow. For example, the large amount of data existing in the literature for displacement flows in uniform channels may be directly extended to the non-uniform channel through scaling of the interface height/speed. Moreover, defining a displacement efficiency parameter,
${\it\eta}$
, as the area fraction behind the front displaced at time
$T$
, we can find
${\it\eta}\approx 1/((1+{\it\alpha}X_{f})V_{f})\approx 1/V_{f}^{\ast }$
. This means that the displacement efficiency is the same for uniform and non-uniform channel flows in the no-back-flow regime. In summary, at time
$T$
we can directly obtain the following relations between the frontal features in uniform and non-uniform channels:
$$\begin{eqnarray}\displaystyle & \displaystyle X_{f}=\frac{-1+\sqrt{1+2{\it\alpha}V_{f}^{\ast }T}}{{\it\alpha}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle V_{f}=\frac{V_{f}^{\ast }}{\sqrt{1+2{\it\alpha}V_{f}^{\ast }T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{f}=h_{f}^{\ast }\sqrt{1+2{\it\alpha}V_{f}^{\ast }T}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\eta}={\it\eta}^{\ast }, & \displaystyle\end{eqnarray}$$
$^{\ast }$
denotes the corresponding value in a uniform channel geometry.
There is an explanation as to why the similarity variables/reductions discussed must exist. The solutions for uniform and non-uniform channel flows are similar in the sense that, as the front advances, different heights of channel are locally achieved and this comes back in some of the similarity solutions. Therefore, provided that an appropriate scaling is adopted, the local base flow velocities that feed into the lubrication model (and therefore the fluxes) will be locally similar. To find the appropriate scaling, for a given
${\it\alpha}X$
, we can use the local channel thickness (
$\hat{D}_{0}(1+{\it\alpha}X)$
) and the local mean flow velocity (
$\hat{V}_{0}/(1+{\it\alpha}X)$
) to rescale the various dimensionless groups, leading to
$$\begin{eqnarray}\displaystyle & \displaystyle h^{\ast }({\it\alpha}X)=h(1+{\it\alpha}X)^{-1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle m^{\ast }({\it\alpha}X)=m(1+{\it\alpha}X)^{2(n_{H}-n_{L})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\chi}^{\ast }({\it\alpha}X)={\it\chi}(1+{\it\alpha}X)^{2n_{H}+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B_{k}^{\ast }({\it\alpha}X)=B_{k}(1+{\it\alpha}X)^{2n_{H}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Re^{\ast }({\it\alpha}X)=Re(1+{\it\alpha}X)^{2(n_{H}-1)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Fr^{\ast }({\it\alpha}X)=Fr(1+{\it\alpha}X)^{-3/2}. & \displaystyle\end{eqnarray}$$
$X$
, confirming that they are similar when the rescaled dimensionless groups (denoted by
$^{\ast }$
in the relations above) are used.
So far, we have been able to define/quantify the characteristics of the no-back-flow regime through a similarity form and by looking at the cases where buoyancy is weak (i.e.
${\it\chi}$
is small). However, our investigation shows that, when
${\it\chi}$
increases, a deviation from the similarity form starts to appear, especially when
${\it\chi}$
exceeds a critical value for which a sustained back-flow regime appears. For example, figure 8 plots the interfaces of figure 6 versus the similarity variable, where it is evident that the interfaces do not collapse for large
${\it\chi}$
when
${\it\alpha}\neq 0$
. Therefore, for strongly buoyant flows in non-uniform channels, it is not possible to extract the behaviour of the interface (e.g. interfacial heights and speeds) directly from the corresponding uniform channel flows. We will quantify the various features of these flows in § 4.4.
Figure 8. The same as figure 6 but plotted versus the similarity variable.
Figure 9. Local base flow velocity profiles and fluxes in a converging channel (
${\it\alpha}=-0.01$
) at
$X=0$
(solid line),
$X=20$
(dashed line) and
$X=50$
(dash-dotted line) using the rescaled dimensionless groups for
${\it\chi}^{\ast }=15$
,
$m^{\ast }=0.1$
. (a) Newtonian displacement: local base flow velocity profiles at for
$h^{\ast }=0.7$
. (b) Newtonian displacement: local fluxes versus
$h^{\ast }=h/(1+{\it\alpha}X)$
. (c) Non-Newtonian displacement: local base flow velocity profiles for
$h^{\ast }=0.7$
,
$B_{H}^{\ast }=2$
,
$B_{L}^{\ast }=5$
,
$n_{H}=n_{L}=0.5$
. (d) Non-Newtonian displacement: local fluxes versus
$h^{\ast }=h/(1+{\it\alpha}X)$
with the parameters of panel (c). The interface slope is zero in all the panels.
Figure 10. Parameter regime in the (
$m,{\it\chi}$
) plane in which there exists multiple fronts of displacements (marked by the shaded area). Elsewhere there is only a single front. In each row,
${\it\alpha}$
is fixed:
$-$
0.01 (a–c), 0.01 (d–f). In each column,
$X_{f}$
is fixed: 5 (a,d), 20 (b,e), 50 (c,f). The insets show the interface at different
$X_{f}$
for the same
${\it\chi}$
and
$m$
.
4.3.1 Single or multiple fronts
Regardless of the displacement flow geometry, we have observed that, for certain parameter sets, the interface has multiple fronts. For a uniform channel Newtonian flow, the competing displacement flow effects are buoyancy (quantified by
${\it\chi}$
), driven by the downhill slope that spreads the interface, and the viscosity ratio (quantified by
$m$
) that sharpens the front. Multiple fronts appear when these physical effects are opposing each another. This is also the case for a non-uniform channel flow, but it may become more interesting as the leading front locally experiences decreasing or increasing buoyancy forces as it advances in a converging or diverging channel, respectively. Through repeated computation of
$\tilde{q}$
for different
$m$
,
${\it\chi}$
,
${\it\alpha}$
and
$X_{f}$
, we can approximate the regime of multiple fronts in the plane of
$(m,{\it\chi})$
for various
${\it\alpha}$
and
$X_{f}$
, as seen in figure 10. To obtain this figure,
$(\partial \tilde{q}/\partial h)(h,{\it\alpha}X)$
(for
$X\in [X_{t},X_{f}]$
) has been very crudely approximated by
$(\partial \tilde{q}/\partial h)(h,{\it\alpha}X_{f})$
, which is then simply inspected for the existence of a second shock. Our simulations confirm that, using this crude approach, we can approximate the boundary between single and multiple fronts. As the leading front progresses forwards (i.e.
$X_{f}$
increases), the shaded region shrinks for a converging channel flow, while it expands for a diverging one. This signifies that it is more (less) likely to find multiple fronts in a diverging (converging) channel as the leading front advances forwards. For example, the insets in figure 10(a–c) for a converging channel flow show that one of the two fronts advancing at short distances disappears, leaving the displacing fluid with a single leading front. Similarly, an interface with a single leading front in a diverging channel flow may split to multiple fronts (although the results are not shown).
4.3.2 Front shape
Regarding the front height, (4.4) is exactly the same equation that would be solved for the hyperbolic conservation law to compute a kinematic shock. However, the front in both uniform and non-uniform channels is not a shock, since diffusive effects of the interface slope are present for
$h\in (0,1+{\it\alpha}X_{t})$
. Finding the approximate shape of the front is straightforward. We should first find
$h_{f}$
from (4.4) and then use
$V_{f}\approx (\partial q/\partial h)(h_{f},{\it\alpha}X_{f})$
. We then shift to a moving frame of reference,
$Z=X-V_{f}T$
, and try to find a steadily travelling solution to (2.22). This relation that we are seeking needs to satisfy
$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}Z}[hV_{f}-q(h,{\it\chi},h_{Z},m,{\it\alpha}X_{f})]=0\quad \Longrightarrow \quad hV_{f}-q(h,{\it\chi},h_{Z},m,{\it\alpha}X_{f})=0. & & \displaystyle\end{eqnarray}$$
Figure 11. Example front shapes in the moving frame of reference, computed from (4.16) for
${\it\chi}=0$
and
$m=1$
at different distances, i.e.
$X_{f}=0$
(solid line),
$X_{f}=20$
(dashed line) and
$X_{f}=50$
(dash-dotted line):
${\it\alpha}=-0.01$
(a), 0 (b), 0.01 (c).
The next step is to simply solve (4.16) numerically for
$h\in (0,h_{f})$
, to find the shape of the front, examples of which are shown in figure 11. This figure shows leading displacement front shapes for fixed
$m=1$
and
${\it\chi}=0$
for various distances (
$X_{f}$
) in both uniform and non-uniform channels. These results clearly show that the frontal region is not a kinematic shock, as it is a region in which the diffusive effects of the interface slope remain significant. Figure 11(a) shows that in a converging channel the frontal diffusive region is shortened as the front advances. Figure 11(b) shows that the frontal region in a uniform channel does not vary with distance as expected. Figure 11(c) for a diverging channel flow illustrates that the transition from the stretched interface region to the frontal region is smoothed out. The front is more diffusive in a diverging channel, i.e. as the front advances the frontal region is extended axially along the channel. Finally, it can be seen in figure 11 that the stretched interface segment also becomes nearly horizontally when
$h$
approaches
$h_{f}$
, which is compatible with the solution of the interface at long times being smooth at the front height.
4.4 Flow regimes for large buoyancy
Regarding the large-scale features of displacement flows, it is critical to identify whether or not a back-flow occurs. The displacement flows without a back-flow are those where the displaced fluid does not travel upstream during the displacement process. Therefore, in order to quantify this, it is important to study the features of the trailing front. We have found that, as the relative buoyancy forces increase, or equally the relative viscous forces decrease, the trailing front presents distinct interesting patterns, which can be easily observed when increasing
${\it\chi}$
from zero. For small
${\it\chi}$
, due to the weak buoyancy, the trailing front does not move upstream during the entire displacement process. In this case, the trailing front seems to be pinned to the upper wall of the uniform/non-uniform channel. We term this flow a no-back-flow regime. As the leading front travels downstream, the thickness of the displacement layer close to the pinned point progressively decreases. Therefore, experimentally, one may expect this point to move at some later time due to molecular diffusion. This effect is obviously absent from our model. For a critical larger value of
${\it\chi}$
, we find a marginal state where the trailing front is at the verge of moving upstream. In this case, the displaced layer has a relatively large finite thickness identified by the trailing front height at
$X\approx 0$
. Meanwhile, the downstream leading front is advected along the channel, leaving behind an apparently stationary residual layer of the displaced layer. This phenomenon has been observed and studied before for uniform pipes and channels, analytically, experimentally and computationally (see e.g. Taghavi et al.
Reference Taghavi, Seon, Wielage-Burchard, Martinez and Frigaard2011). Curiously, the value of the critical
${\it\chi}={\it\chi}_{c}$
at which the stationary interface flow regime occurs is independent of
${\it\alpha}$
. For even larger values of
${\it\chi}$
, the trailing front starts to move backwards. For uniform and converging channels, we have observed that the trailing front advances continuously upstream, demonstrating a sustained back-flow regime. However, for any given large
${\it\chi}$
in a diverging channel, although the trailing front initially starts to move upstream, it slows down and eventually stops at some location. We have christened this flow an eventually stationary interface flow regime. When the trailing front propagates backwards in a diverging channel, it is in fact moving in a progressively narrower channel, for which the local buoyancy stresses are progressively smaller. The back-flow is able to resist the imposed flow until it reaches a critical upstream location where the viscous stresses perfectly balance the local buoyancy stresses. Thus, the back-flow stops. After this critical distance is reached, the flow behaviours become similar to that of a stationary interface flow regime.
The stationary interface flow identifies the transition between the no-back-flow and the sustained back-flow regimes for uniform and converging channels. It equally identifies the transition between the no-back-flow and the eventually stationary interface flow regimes for diverging channels. The existence of a displaced fluid layer that is apparently stationary with a constant flow rate of displacing heavy fluid implies that the interfacial speed is zero and the flux is unity at the interface. At long times, the interface is expected to elongate so that the effect of the slope of the interface in all regions may be negligible except local to the advancing frontal regions. Therefore, to find the condition of the stationary interface, we can rely on the approximate hyperbolic part of (2.22) by setting
$h_{X}=0$
. Taghavi et al. (Reference Taghavi, Seon, Wielage-Burchard, Martinez and Frigaard2011, Reference Taghavi, Alba and Frigaard2012a
) and Moyers-Gonzalez et al. (Reference Moyers-Gonzalez, Alba, Taghavi and Frigaard2013) have quantified the conditions for the occurrence of stationary flows for various parameters, e.g. the viscosity ratio, the buoyancy number and even viscoplastic rheology, for uniform channels and pipes. Our investigation shows that the conditions for the occurrence of stationary interface flow regimes is not affected when the flow geometry becomes slightly non-uniform.
Let us consider the case of displacement flows in a diverging channel separately. For a given viscosity ratio,
$m$
, and any given
${\it\chi}>{\it\chi}_{c}$
, our simulations show that there is a critical upstream distance where the back-flow stops moving, demonstrating an eventually stationary interface flow regime. When the back-flow stops, the interface is elongated, varying in the interval
$h\in [0,1+{\it\alpha}X_{es}]$
, where the subscript es refers to eventually stationary. The interface at
$X_{es}$
has a finite height denoted by
$h_{es}$
. Two conditions must be simultaneously met for the interface to eventually become stationary, i.e.
$$\begin{eqnarray}\displaystyle \int _{0}^{h_{es}}\frac{\partial \tilde{q}(h,{\it\alpha}X)}{\partial h}\,\text{d}h=\tilde{q}(h_{es},{\it\alpha}X_{es})=1, & & \displaystyle\end{eqnarray}$$
and the approximate relation for the interface speed at
$h_{es}$
,
$$\begin{eqnarray}\displaystyle V_{es}\approx \frac{\partial \tilde{q}}{\partial h}(h_{es},{\it\alpha}X_{es})=0. & & \displaystyle\end{eqnarray}$$
Equations (4.17) and (4.18) can be solved to obtain
$X_{es}$
for given
$m$
and
${\it\chi}$
. Although in general we need to numerically find the conditions for the eventually stationary interface flow regime, for
$m=1$
we can derive an analytical relationship between a given
${\it\chi}_{es}$
and the critical upstream distance,
$X_{es}$
, where the back-flow stops:
$$\begin{eqnarray}\displaystyle {\it\chi}_{es}=12\frac{\sqrt{2}(1+\sqrt{2})}{(-2+\sqrt{2})(1+{\it\alpha}X_{es})^{3}},\quad {\it\alpha}>0\text{ and }X_{es}<0. & & \displaystyle\end{eqnarray}$$
After some algebra, for large viscosity ratios
$m\gg 1$
, we can also find
$$\begin{eqnarray}\displaystyle {\it\chi}_{es} & = & \displaystyle -3\frac{\sqrt{5}+1}{(-2+\sqrt{5})(1+{\it\alpha}X_{es})^{3}(\sqrt{5}-3)^{3}}-9\frac{1}{(1+{\it\alpha}X_{es})^{3}(\sqrt{5}-3)^{4}m}+O(m^{-2}),\qquad \nonumber\\ \displaystyle & & \displaystyle \quad {\it\alpha}>0\text{ and }X_{es}<0.\end{eqnarray}$$
Finally, for small viscosity ratios
$m\ll 1$
, we arrive at
$$\begin{eqnarray}\displaystyle {\it\chi}_{es}=\frac{1}{6}\frac{(36^{2/3}+102\sqrt[3]{6}m^{2/3}+63m^{1/3})\sqrt[3]{6}}{({\it\alpha}X_{es}+1)^{3}}+O(m),\quad {\it\alpha}>0\text{ and }X_{es}<0. & & \displaystyle\end{eqnarray}$$
From the equations above, it is clear that
${\it\chi}_{es}\rightarrow {\it\chi}_{c}$
as
$X_{es}\rightarrow 0$
. Thus, we can simplify and write
$$\begin{eqnarray}\displaystyle X_{es}\approx -\!\left(1-\left(\frac{{\it\chi}_{c}}{{\it\chi}_{es}}\right)^{1/3}\right){\it\alpha}^{-1},\quad {\it\alpha}>0\text{ and }{\it\chi}_{es}>{\it\chi}_{c}, & & \displaystyle\end{eqnarray}$$
where
${\it\chi}_{c}$
is a constant that can be found as
${\it\chi}_{c}=3.3$
for
$m\ll 1$
,
${\it\chi}_{c}=69.9$
for
$m=1$
, and
${\it\chi}_{c}=92.2$
for
$m\gg 1$
. First of all, (4.22) quantifies that, for any large value of
${\it\chi}_{es}$
, there exists a critical upstream distance where the interface eventually becomes stationary. Second, this equation shows that, for very large values of
${\it\chi}_{es}$
, the stopping distance of the trailing front varies like
$X_{es}\sim {\it\alpha}^{-1}$
.
Figure 12. (a) Variation of
${\it\chi}_{es}$
versus
$X_{es}$
for large viscosity ratios (i.e.
$m=10$
) showing the numerical solution of (4.17) and (4.18) (solid line) and (4.20) (dashed line). (b) Variation of
${\it\chi}_{c}$
versus
$X_{es}$
for small viscosity ratios (i.e.
$m=0.1$
) showing the numerical solution of (4.17) and (4.18) (solid line) and (4.21) (dashed line). (c) Variation of
${\it\chi}_{c}$
versus
$X_{es}$
for large viscosity ratios (i.e.
$m=1$
) using (4.19). Black circles show the simulation results. The insets are explained in the text.
Figure 12 shows the variation of
${\it\chi}_{es}$
versus
$X_{es}$
for various viscosity ratios. Figure 12(a) shows that the prediction of analytical equation (4.20) and the solution of (4.17) and (4.18) are in perfect agreement with each other, and also with the simulation data (black circles) for large
$m$
. Figure 12(b) shows that the prediction of analytical equation (4.21) provides a good approximation to the solution of (4.17) and (4.18) matching with the simulation data for small
$m$
. Figure 12(c) shows that (4.19) matches the simulation data for
$m=1$
. The inset in figure 12(a) depicts the trailing front evolution from the simulations, showing that the trailing front slows down until it stops. The inset in figure 12(b) plots the variation of the trailing front velocity (from the simulations) eventually reaching zero for the data point marked by the arrow. Note that, to obtain the simulation data, for certain cases it takes a very long time (e.g.
$T\sim O(10^{3})$
) for the back-flow to stop completely. Within this time, for a miscible displacement, molecular diffusion may come into play and modify the flow regime, the characterization of which is beyond the scope of this work.
5 Non-Newtonian fluids
Before we proceed to provide the results for non-Newtonian displacements, two points are worth mentioning, as follows.
-
(i) Overall, the regime classifications discussed for Newtonian displacements are also intuitively valid for non-Newtonian fluids. For small
${\it\chi}$
, we have the no-back-flow regime; and for sufficiently large
${\it\chi}$
, we pass through the stationary interface flow regime and enter the sustained back-flow regime (for uniform and converging channels) or eventually stationary interface flow regime (for diverging channels). Similar to Newtonian displacements, various flow features, such as the condition for the existence of multiple fronts, the conditions for the stationary and the eventually stationary interface flow regimes, front shapes, front heights and speeds, etc., can be obtained for shear-thinning and yield-stress fluid displacements. We will not present these results, since we will be focusing on presenting other interesting effects relevant to non-Newtonian features. -
(ii) The superposition of the interface heights for the no-back-flow regime can be obtained for shear-thinning fluids (although the results are not shown), but for yield-stress fluids it is not possible, even in the case of uniform channels. Therefore, one has to perform simulations and cannot simply obtain the behaviour of the interface at long times through a similarity form.
5.1 Shear-thinning effects
The Herschel–Bulkley model considers shear-thinning and yield-stress features of non-Newtonian fluids. Here, we start by exploring only shear-thinning effects,
$B_{k}=0$
. Many features of the shear-thinning fluid displacement may be predictable, even for non-uniform geometry. For example, for both uniform and non-uniform channels, we observe that, as
$n_{L}$
(
$n_{H}$
) decreases, the front height increases (decreases) for fixed
$n_{H}$
(
$n_{H}$
). This is similar to effects of the viscosity ratio parameter, as varying the power-law indices makes one fluid progressively more or less viscous. To observe a significant effect of the geometry on the flow, we concentrate on a converging channel flow with no buoyancy. Figure 13 shows the interface evolution with time, showing a shear-thinning fluid displacing a Newtonian fluid at
${\it\chi}=0$
,
$m=1$
and three power-law indices, i.e.
$n_{H}=0.75,0.5,0.25$
. Figure 13 is interesting: as the channel narrows down, the leading front is naturally advancing faster and therefore it experiences a larger shear rate. Consequently, the displacing fluid becomes progressively less viscous so that the relative local ratio of the viscosity will drop further, causing the leading front to slump more and move even faster. This effect becomes magnified especially at longer distances: as can be seen in figure 13, the front is moving very fast for the lowest value of
$n_{H}$
at large distances.
5.2 Yield-stress effects
In this section, we discuss the general effects of a yield stress present in the displaced fluid on the displacement flow. Some of features might be common between uniform and non-uniform channel flows. For example, the displacement efficiency is usually enhanced through the effects that make the displacing fluid more viscous. For instance, it may be expected that increasing
$B_{H}$
improves the efficiency of the displacement, thanks to the enhanced effective viscosity, and increasing
$B_{L}$
makes the displacement less efficient. Another physical phenomenon that is expected is that, for larger
$B_{L}$
, it is possible to have a static residual wall layer, in which fluid
$L$
occupying the upper part of the channel is not displaced. Static residual wall layers have been observed and studied in uniform geometries, for a static situation as well as a transient displacement flow (see e.g. Allouche et al.
Reference Allouche, Frigaard and Sona2000; Frigaard, Leimgruber & Scherzer Reference Frigaard, Leimgruber and Scherzer2003). We will discuss static residual wall layer solutions further in § 5.2.1.
Some effects are definitely harder to predict. For example, varying
$m$
,
${\it\chi}$
and
$B_{k}$
at the same time may result in behaviours that might not be expected beforehand. Some other features, such as the existence of multiple fronts, the onset of the back-flow regime and the critical stopping distance of the back-flow, may not be known beforehand. Our model provides a tool to approximate these features for two generalized Newtonian fluids; however, employing our model to derive all these conditions is lengthy and beyond the scope of this paper. In addition, it should be mentioned that, regardless of the flow geometry, there seems to be two critical
${\it\chi}$
for a yield-stress fluid: first, for
${\it\chi}={\it\chi}_{c}$
the displaced fluid is about to move backwards; second, for
${\it\chi}={\it\chi}_{c}^{\prime }$
the static residual layers of the displaced fluid are about to be displaced forwards (i.e. when the absolute value of the upper wall shear stress is equal to the yield stress). It is interesting that for
${\it\chi}_{c}^{\prime }<{\it\chi}<{\it\chi}_{c}$
the lighter displaced fluid is not displaced. We leave the in-depth study of these critical
${\it\chi}$
for a future work.
The most interesting case to look at would be perhaps the effect of non-uniform geometry on the displacement of a viscoplastic fluid by a Newtonian fluid. We also set
$n_{L}=1$
, for simplicity, so that the displaced fluid is a Bingham fluid. Owing to the yield stress, these fluids are shear-thinning while no additional power-law behaviour is imposed.
Figure 13. Examples showing a shear-thinning fluid displacing a Newtonian fluid for
${\it\chi}=0$
,
$m=1$
and
$T=0,3,\ldots ,30$
for
${\it\alpha}=-0.01$
:
$n_{H}=0.75$
(a), 0.5 (b), 0.25 (c).
Figure 14. Interface evolution in time for
$T=0,3,\ldots ,30$
with
$m=1$
,
${\it\chi}=0$
and
$n_{L}=1$
. In each row,
${\it\alpha}$
is fixed: 0 (a–c),
$-$
0.01 (d–f), 0.01 (g–i). In each column,
$B_{L}$
is fixed: 5 (a,d,g), 20 (b,e,h), 40 (c,f,i). The displacing fluid is Newtonian.
Figure 14 shows examples of Bingham displacements for a fixed viscosity ratio (
$m=1$
) and no buoyancy (
${\it\chi}=0$
). In each row of figure 14 the yield stress of the displaced fluid increases as
$B_{L}=5$
(a,d,g), 20 (b,e,h) and 40 (c,f,i). The heights of the maximal static residual wall layers, approximated by the theory of § 5.2.1, are superimposed as dash-dotted lines on all the plots, showing reasonable agreement. Figure 14(a–c) shows the results for a uniform channel flow. In (d–f), for a converging channel flow, we can see that at
$B_{L}=5$
there is no static residual wall layer (figure 14
d). By increasing the yield stress to
$B_{L}=10$
there exists a relatively thick static residual wall layer, at small distances (figure 14
e). We observe that the predictions of the static residual wall layer thickness at smaller
$X$
perfectly matches the simulation results. We also note that the height of the leading front decreases as the displacement flow develops; however, the height of the channel decreases at a faster rate, so that static residual wall layer thickness has to decrease. This decline continues until a critical distance at which the static residual wall layer completely disappears. We quantify this distance in § 5.2.1. Figure 14(f) shows more or less the same behaviour as observed in figure 14(e) but for an initially thicker static residual wall layer. Figure 14(g–i) displays the results for a diverging channel flow. Figure 14(g) is noteworthy in that there is initially no static residual wall layer but at a certain distance, as the dash-dotted line shows, a progressively growing static residual wall layer appears. In figure 14(h,i), we observe that, starting around
$T\sim O(1)$
, there exists a static residual wall layer, the thickness of which is progressively increasing.
Figure 15. Examples of Bingham fluids (
$B_{L}=20$
,
$n_{L}=1$
) displaced by Newtonian fluids for
${\it\chi}=200$
,
$m=1$
and
$T=0,3,\ldots ,30$
;
${\it\alpha}=0$
(a),
$-$
0.01 (b), 0.01 (c).
Before we proceed further, we can draw qualitative comparisons between our results and those available in the literature with regard to static layers in non-uniform channels, e.g. with the results of Roustaei & Frigaard (Reference Roustaei and Frigaard2013) (see also Roustaei & Frigaard Reference Roustaei and Frigaard2015; Roustaei, Gosselin & Frigaard Reference Roustaei, Gosselin and Frigaard2015), who studied the onset and characteristics of static layers through extensive Stokes flow computations of a single Bingham fluid in non-uniform channels. (i) Our study and theirs both find that the combination of rheology and channel non-uniformity causes the formation (or destruction) of static layers (see figure 14
f,e and further discussions in § 5.2.1). Also, both works find that static layers intuitively appear first at the widest part of the channel. (ii) Roustaei & Frigaard (Reference Roustaei and Frigaard2013) discovered that static layers are found at a certain depth of the wall perturbation. Comparably, we find that static layers appear at a certain length of a diverging channel, since, as the front advances, larger depths of the channel are achieved (see e.g. figure 14
f). (iii) They showed that, when the static layers are formed, their thickness grows as the depth of the wall perturbation increases. We also find that, as
$X$
increases, the static layer thickness increases. (iv) Roustaei & Frigaard (Reference Roustaei and Frigaard2013) found that the product
$\hbar {\it\delta}$
(with
${\it\delta}$
the aspect ratio and
$\hbar$
the maximal depth of the wall perturbation) is the relevant geometrical parameter that governs the onset of static layers. For example, for sufficiently large
$\hbar {\it\delta}$
(at fixed Bingham number), there exist static layers. Considering the lubrication rescaling, our comparable geometrical parameter is
${\it\alpha}X$
, which has the same effect as
$\hbar {\it\delta}$
. (v) Finally, Roustaei & Frigaard (Reference Roustaei and Frigaard2013) provided relationships of
$\hbar {\it\delta}$
as a function of the Bingham number to portray the existence of static layers. We can also quantify the existence of these layers for our buoyant two-fluid system, examples of which are relations (5.12) and (5.13), discussed in § 5.2.1.
One may wonder about the effect of combining a yield stress of the displaced fluid with buoyancy. Our investigation has shown that, for small buoyant forces (i.e. small
${\it\chi}$
), the buoyancy does not modify the flow features. However, for a relatively large
${\it\chi}$
, the flow can present stimulating patterns. For example, figure 15 shows examples of Bingham displacements (
$B_{L}=20$
) for
${\it\chi}=200$
and
$m=1$
, in the three channel geometries. Figure 15(a) shows that there exists a static residual wall layer and there are also two fronts travelling with different constant speeds. In figure 15(b), there is initially a static residual wall layer at small distances but it eventually disappears. There are also two fronts: the faster one initially seems to have more or less a constant height and a constant speed; however, at some distances it suddenly starts to feel the presence of the converging upper wall. At this point the front height suddenly drops and the front speed increases significantly. The second front follows the first one with progressively smaller speed; consequently, the second slows down and eventually stops moving. Figure 15(c) shows that there exists a static residual wall layer in a diverging channel and there are initially two fronts. The lower front (i.e. closer to the lower wall) is initially moving faster. Later, the second front moves faster and merges with the first one to form a single front. Figure 15 also reveals that the theory of § 5.2.1 generally overestimates the thickness of the static layer for large
${\it\chi}$
.
5.2.1 Static residual wall layers
Static residual wall layers have been first analysed in the displacement flow context by Allouche et al. (Reference Allouche, Frigaard and Sona2000), who considered displacement flows in a vertical uniform channel. In their setting, the interface slope does not enter the equations for the interface propagation and the lubrication displacement model is hyperbolic. For parameters that admit a static residual wall layer, Allouche et al. (Reference Allouche, Frigaard and Sona2000) found that the interface approaches the maximal static layer thickness at long times. In our model, away from the moving fronts, the flow at long times is governed by the hyperbolic model by abandoning the slope of the interface. Therefore, we expect that the interface will also approach the maximal layer thickness at long times. Here, we outline the calculation of the maximal layer thickness in non-uniform channel geometry. We rely on a similar procedure to that presented in Allouche et al. (Reference Allouche, Frigaard and Sona2000) and Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009).
First, let us suppose that the light displaced fluid is static for
$y\in [h,1+{\it\alpha}X]$
, for some
$h\in (0,1+{\it\alpha}X]$
. The flow in the heavy fluid layer is thus governed by
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial y}{\it\tau}_{H,Xy}=\frac{\partial P_{0}}{\partial X},\quad y\in (0,h), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u(0)=u(h)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \int _{0}^{h}u(y)\,\text{d}y=1. & & \displaystyle\end{eqnarray}$$
On integrating to find the velocity field and then integrating across the heavy fluid layer, (5.3) becomes
$$\begin{eqnarray}\displaystyle \frac{1}{2}\left(\frac{\tilde{B}_{H}}{{\it\zeta}}\right)^{1/n_{H}}h^{2}\left[\frac{(1-{\it\zeta})^{1+1/n_{H}}}{1+\displaystyle \frac{1}{n_{H}}}-\frac{(1-{\it\zeta})^{2+1/n_{H}}}{\left(1+\displaystyle \frac{1}{n_{H}}\right)\left(2+\displaystyle \frac{1}{n_{H}}\right)}\right]=1, & & \displaystyle\end{eqnarray}$$
where
$\tilde{B}_{H}=B_{H}/{\it\kappa}_{H}$
, and
${\it\zeta}$
is related to
$\partial P_{0}/\partial X$
by
$$\begin{eqnarray}\displaystyle \frac{\partial P_{0}}{\partial X}=-\frac{2B_{H}}{{\it\zeta}h}. & & \displaystyle\end{eqnarray}$$
It can be shown that (5.4) has a unique root that lies between
$0$
and
$1$
, which can be found numerically. Thus, we may write
${\it\zeta}={\it\zeta}(h,\tilde{B}_{H},n_{H})$
. We suppose that
${\it\chi}$
is sufficiently small that the stress within the light fluid layer at the upper wall has the same sign as the interfacial stress. The wall shear stress at the upper wall is
$-(1+{\it\alpha}X-0.5h)(-\partial P_{0}/\partial X)+{\it\chi}(1+{\it\alpha}X-h)$
and the initial assumption of a static layer is only valid provided that
$$\begin{eqnarray}\displaystyle \left|\left(1+{\it\alpha}X-\frac{h}{2}\right)\left(-\frac{\partial P_{0}}{\partial X}\right)-{\it\chi}(1+{\it\alpha}X-h)\right|\leqslant B_{L}. & & \displaystyle\end{eqnarray}$$
We denote the minimal value of
$h$
for which (5.6) holds by
$h_{min}$
to define a maximal static layer by
$$\begin{eqnarray}\displaystyle Y_{static}=1+{\it\alpha}X-h_{min}. & & \displaystyle\end{eqnarray}$$
Substituting
$\partial P_{0}/\partial X$
from (5.5) into (5.6), we find
$$\begin{eqnarray}\displaystyle \frac{(2+2{\it\alpha}X-h_{min})\tilde{B}_{Y}}{{\it\zeta}(h_{min},\tilde{B}_{H},n_{H})h_{min}}-(1+{\it\alpha}X-h_{min})\tilde{{\it\chi}}=1. & & \displaystyle\end{eqnarray}$$
Therefore, the computation of the maximal residual wall layer yields a solution that depends only on five parameters, i.e.
$n_{H}$
,
$\tilde{B}_{H}=B_{H}/k_{H}$
,
$\tilde{B}_{Y}=B_{H}/B_{L}$
,
$\tilde{{\it\chi}}={\it\chi}/B_{L}$
and
${\it\alpha}X$
.
The critical condition for which any static residual wall layer exists occurs when
$h_{min}\rightarrow 1+{\it\alpha}X$
, and it is seen to be completely independent of the buoyancy ratio
$\tilde{{\it\chi}}$
. This is given by
$$\begin{eqnarray}\displaystyle \frac{\tilde{B}_{Y}}{{\it\zeta}(1+{\it\alpha}X,\tilde{B}_{H},n_{H})}=1. & & \displaystyle\end{eqnarray}$$
Figure 16. Maximal static residual wall layer thickness,
$Y_{static}$
, with contours spaced at intervals
$Y_{static}=0.05$
for
$\tilde{{\it\chi}}=1$
and
$n_{H}=1$
. In each row,
${\it\alpha}$
is fixed:
$-0.01$
(a–c), 0.01 (d–f). In each column,
$X$
is fixed: 10 (a,d), 50 (b,e), 90 (c,f).
The maximum static residual wall layer
$Y_{static}$
varies with the parameters
$\tilde{B}_{Y}$
and
$\tilde{B}_{H}$
, as seen in figure 16 for converging and diverging channels at three distances. Using the shaded area, we have marked on this figure the limit where static residual wall layers cannot exist. Figure 16(a–c) shows that the region of no static residual wall layer expands as
$X$
increases, i.e. when the displacing fluid progresses forwards. This means that it is more likely for the static residual wall layers to disappear at large distances in a converging channel. Figure 16(d–f) shows that in a diverging channel the region of no static residual wall layer slightly shrinks as
$X$
increases. This implies that it is more likely for the static residual wall layers to appear at large distances in a diverging channels.
Note that the limit
$B_{H}\rightarrow 0$
is indeterminate in the above calculation, and so must be treated separately. The relation between flow rate and pressure drop is easily found for any static
$h$
, which we insert into (5.6), furnishing the minimal
$h$
(or maximal static layer):
$$\begin{eqnarray}\displaystyle \frac{2(1+{\it\alpha}X)-h_{min}}{h_{min}^{2n_{H}+1}}\left(\frac{4n_{H}+2}{n_{H}}\right)^{n_{H}}-(1+{\it\alpha}X-h_{min})\tilde{\tilde{{\it\chi}}}-\tilde{B}_{L}=0, & & \displaystyle\end{eqnarray}$$
where
$\tilde{\tilde{{\it\chi}}}={\it\chi}/{\it\kappa}_{H}$
and
$\tilde{\tilde{B}}_{L}=B_{L}/{\it\kappa}_{H}$
. The above equation contains interesting findings. For example, one can postulate that, in a diverging channel at
$X\rightarrow \infty$
,
$h_{min}$
may reach a constant value, which we can calculate as
$$\begin{eqnarray}\displaystyle h_{min}\approx \left(\frac{2}{\tilde{\tilde{{\it\chi}}}}\left(\frac{4n_{H}+2}{n_{H}}\right)^{n_{H}}\right)^{1/(2n_{H}+1)}. & & \displaystyle\end{eqnarray}$$
Therefore, for any
$\tilde{\tilde{{\it\chi}}}$
, there exists an
$h_{min}$
, meaning that there exists a static residual wall layer at least at large distances. Equation (5.10) also provides other results. For example, it can be shown that, in a diverging channel (i.e.
${\it\alpha}>0$
), even for small values of
$\tilde{\tilde{B}}_{L}$
, there is a distance at which the static residual wall layer starts to appear:
$$\begin{eqnarray}\displaystyle X_{static}^{appear}=\frac{1}{{\it\alpha}}\left(\sqrt{\frac{1}{\tilde{\tilde{B}}_{L}}\left(\frac{4n_{H}+2}{n_{H}}\right)^{n_{H}}}-1\right),\quad 0<\tilde{\tilde{B}}_{L}<\left(\frac{4n_{H}+2}{n_{H}}\right)^{n_{H}}. & & \displaystyle\end{eqnarray}$$
On the other hand, in a converging channel (i.e.
${\it\alpha}<0$
), even for large values of
$\tilde{\tilde{B}}_{L}$
, there is a distance at which the static residual wall layer has to disappear:
$$\begin{eqnarray}\displaystyle X_{static}^{disappear}=\frac{1}{{\it\alpha}}\left(\sqrt{\frac{1}{\tilde{\tilde{B}}_{L}}\left(\frac{4n_{H}+2}{n_{H}}\right)^{n_{H}}}-1\right),\quad \tilde{\tilde{B}}_{L}>\left(\frac{4n_{H}+2}{n_{H}}\right)^{n_{H}}. & & \displaystyle\end{eqnarray}$$
We note that
$X_{static}^{disappear}\rightarrow {\it\alpha}^{-1}$
as
$\tilde{\tilde{B}}_{L}\rightarrow \infty$
, which is consistent with the assumptions of our model.
These aforementioned findings are in agreement with the simulation results in figure 14(e–g).
6 Discussion and conclusions
We have considered the viscous limit of a miscible, significantly buoyant displacement flow in non-uniform channel geometry. We have used a lubrication/thin-film approximation to study fluids of generalized Newtonian type. We can summarize our findings as follows.
For small buoyancy, we find a no-back-flow regime. We have analysed principally the long-time behaviours, where the displacement flow is characterized by two intervals. In the first interval, the trailing front of the interface is pinned to the upper wall as the interface gradually stretches. The lower part of the interface forms a leading front that is sharp and that moves downstream. The leading front propagates in time with a constant speed in uniform channels, an increasing speed in converging channels, and a decreasing speed in diverging channels. Correspondingly, the height of the interface at the front is constant for uniform channel flows, typically decreasing for a converging channel flow, and usually increasing for a diverging channel flow. By consideration of the associated hyperbolic problem, we are able to directly determine the front heights and front velocities for this regime. Quite interesting is that the solutions of the displacement flow at long times in all the geometries considered converge to the similarity form
$(X+{\it\alpha}X^{2}/2)/T$
.
We have also quantified other characteristics of the no-back-flow regime. It is possible for the displacing fluid to exhibit more than one leading front. We quantify this feature through a regime map versus the buoyancy number and the viscosity ratio for converging and diverging channels. As the displacement flow develops, it is progressively less (more) likely to have multiple fronts in a converging (diverging) channel. Another feature is the presence of diffusive effects of the interface slope at the front, which do not vary with distance for a uniform channel, decrease with distance in a converging channel, and expand with distance in a diverging channel.
The appearance of back-flow is a characteristic of the flows with large buoyancy, independent of the flow geometry. For these flows, we have classified the flows using the stationary interface flow regime, the sustained back-flow regime (only for uniform and converging channels), and the eventually stationary interface flow regime (only for diverging channels). The latter is remarkable, since we find that the displaced fluid moves backwards but stops at a critical distance with respect to the initial interface. We have provided predictions of this critical distance by consideration of the associated hyperbolic problem, agreeing well with the simulation results. Nevertheless, the approximate equation from the associated hyperbolic problem cannot provide accurate predictions of the leading front heights of these flows.
Regarding the shear-thinning features, the case of our interest was a shear-thinning fluid displacing a Newtonian fluid in a converging channel. When the power-law index is low, the displacing fluid undergoes a large shear rate at larger distances and it thins relatively more rapidly, resulting in a significant increase in the front speed and a significant decrease in the front height.
When a yield stress is present in the displaced fluid, it is possible for there to exist completely static residual wall layers. Our theoretical analysis delivers the maximal static residual wall layer thickness for displacement flows, which may be significantly affected by the non-uniformity of the geometry. In particular, in a converging channel and for a given parameter set, it is possible for an initially significantly thick static residual wall layer to disappear at a critical distance. In diverging channels, we have quantified that, for any non-zero value of the yield stress, the static residual wall layer starts to appear at a critical distance. In addition, a combination of buoyancy and a large yield stress in the displaced fluid may result in the appearance of a second front moving downstream. This front may develop at a different speed from that of the first leading front in a uniform channel, it may merge with the first leading front at some distance in a diverging channel, or it may stop moving in a converging channel. We recall that our analysis of the static residual wall layer must be seen in the light of providing mainly qualitative understanding rather than definitive answers to the question of where the boundaries of static and yielded regions are.
It is appropriate to discuss some of the important differences between the current problem and the parallel-wall, channel problem studied in Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009). First of all, while developing a classical lubrication model is the main approach for the two works, the current study extends the previous one by including another key dimensionless parameter that is geometrical (i.e.
${\it\alpha}X$
). The analysis is directed towards understanding where this parameter (for small
${\it\alpha}$
) may significantly affect displacement flows. For small buoyancy, a similarity solution form is expectedly generalized and extended to include the geometrical parameter. In a sense, the displacement fronts in these flows simply advance in a progressively narrower/wider channel with no significant or unexpected effect of the geometrical parameter. In contrast, at large buoyancy, while no similarity solution is found, the geometrical parameter becomes more relevant and results in interesting findings (e.g. controlling the trailing front position in a diverging channel). The reason is that, for these flows, the front locally experiences different buoyant forces depending on the spatial location and this modifies the interface behaviour globally. Lastly, the full parabolic problem must be computationally analysed to extract interface behaviours in non-parallel channels, whereas the associated hyperbolic problem is sufficient to directly deliver long-time flow features in parallel channels.
Finally, our model and its results can be validated with experiments and CFD computations. Experiments will be straightforward, although aspect ratios and buoyant effects should be large at the same time, which may impose restrictions on experimental apparatus dimensions. Computations may exhibit more challenges that must be overcome. Note that our model is focused on long-time behaviours (implying large aspect ratios) and include complex fluid rheologies, increasing the number of dimensionless parameters. While three-dimensional computations of miscible displacement flows of Newtonian fluids are common on large parallel machines (e.g. John et al. Reference John, Oliveira, Heussler and Meiburg2013; Hallez & Magnaudet Reference Hallez and Magnaudet2009), existing computational speeds restrict these studies to low/moderate aspect ratios and impose limits on the ranges of dimensionless parameters. Despite the challenges, similar displacement flow models (for uniform geometries) have been numerically and experimentally validated (see e.g. Taghavi et al. Reference Taghavi, Alba, Seon, Wielage-Burchard, Martinez and Frigaard2012c ).
There are possible interesting extensions of the current work, which we leave for future consideration, e.g. including surface tension forces in the model and studying an immiscible version of the displacement flow studied, or studying fluid flow instabilities and other fluid phenomena, which may result in the breakdown of the underlying model assumptions.
Acknowledgements
This research has been carried out at Université Laval, supported financially by the NSERC Discovery Grant. R.M. also acknowledges additional financial support for his postdoctoral programme from the Aluminium Research Centre (REGAL), Université Laval. This support is gratefully acknowledged. The research has been enabled in part by support provided by Calcul Québec. We wish to acknowledge useful discussions with I. A. Frigaard and K. Alba. We also thank the reviewers for their constructive comments.
Appendix A. Alternative scale for viscosity
Based on the scaling balanced between the two fluids, an alternative scale for viscosity is presented here, which may facilitate the understanding of increasing/ decreasing viscosity in the displacing/displaced fluid. While a mean density has already been defined (i.e.
$\hat{{\it\rho}}=(\hat{{\it\rho}}_{H}+\hat{{\it\rho}}_{L})/2$
), we define a mean viscosity for which a choice would be
$$\begin{eqnarray}\displaystyle \hat{{\it\mu}}=\sqrt{\hat{{\it\mu}}_{H}\hat{{\it\mu}}_{L}}. & & \displaystyle\end{eqnarray}$$
Therefore, we can make the Navier–Stokes equations dimensionless, the same as before, except that we should use
$\hat{{\it\mu}}\hat{V}_{0}/\hat{D}_{0}$
to scale pressure and stresses. The mean viscosity is then used to define the Reynolds number in (2.3). Using the new alternative scaling, the ratios of viscosity in the heavy and light layers are
$\hat{{\it\mu}}_{H}/\hat{{\it\mu}}=1/\sqrt{m}$
and
$\hat{{\it\mu}}_{L}/\hat{{\it\mu}}=\sqrt{m}$
, respectively, where
$m$
is the viscosity ratio parameter in the main text. Similarly, the new Bingham numbers are defined as
$$\begin{eqnarray}\displaystyle B_{k}=\frac{\hat{{\it\tau}}_{k,Y}}{\sqrt{\hat{{\it\kappa}}_{H}\hat{{\it\kappa}}_{L}}\left[\displaystyle \frac{\hat{V}_{0}}{\hat{D}_{0}}\right]^{(n_{L}+n_{H})/2}}. & & \displaystyle\end{eqnarray}$$
