1. Introduction
Complete removal of drilling fluids and sealing the space between rock formations and steel casing with a cement slurry is a critical step in the construction of water, geothermal, hydrocarbon and waste storage wells (Guillot et al. Reference Guillot, Hendriks, Callet and Vidick1990). This process is called primary cementing and its foremost purpose is to provide zonal isolation (Nelson & Guillot Reference Nelson and Guillot2006; Lavrov & Torsæter Reference Lavrov and Torsæter2016), which means that the different geological strata penetrated by the well are unable to communicate hydraulically along the completed well. This objective is not always achieved. For example,
$28.5\,\%$
of wells drilled in British Columbia between 2010 and 2018 (Trudel et al. Reference Trudel, Bizhani, Zare and Frigaard2019) and
$28{-}32\,\%$
of wells drilled in Southeastern Saskatchewan (Baillie et al. Reference Baillie, Risk, Atherton, O’Connell, Fougère, Bourlon and MacKay2019), were reported to emit gas. The leakage rates in such surveys are generally heavily skewed towards low rates (
${\lessapprox}1$
m
$^3$
day–1), but nevertheless indicate that the well integrity is often compromised. Estimates of leakage in other regions globally may be significantly lower. This can be because the wells considered are not primarily gas wells. In some jurisdictions data are either not publicly accessible or not recorded.
The primary cementing operation comprises two main steps; see figure 1. After drilling the borehole and removing the drill pipe, a steel casing is inserted into the drilled well that is filled with drilling fluid. Then, a sequence of fluids are pumped from the surface through the casing, travelling down to the bottom of the well and then returning upwards in the annular space between the wellbore wall and the outer surface of the casing. To reduce any chemical incompatibility between the drilling fluid and the cement slurry, preflushes (called washes and spacers) are commonly pumped ahead of the cement slurry (Nelson & Guillot Reference Nelson and Guillot2006; Guillot, Desroches & Frigaard Reference Guillot, Desroches and Frigaard2007). Therefore, the process relies upon a sequence of fluid--fluid displacement flows to be effective.
Schematic of primary cementing operations: (a) inserting the casing into the drilled well, (b) injecting the cement slurry into the annulus.
The flow challenges are significant, as the displacement takes place along many 100’s of metres of a narrow eccentric annulus of changing orientation. The drilling mud has yield stress and is shear thinning, as with many spacer fluids and the cement slurry. Large density differences, together with changing orientation of the well mean that buoyancy is a significant effect, despite a typically small annular gap (
$0.5{-}5$
cm). The yield stress of the drilling mud is essential as part of the drilling process, but can cause problems later. It may resist displacement, allowing the mud to bridge between inner and outer walls of the annulus (a mud channel), remaining after cementing (McLean, Manry & Whitaker Reference McLean, Manry and Whitaker1967). Less severe, the yield stress may not be exceeded at the walls allowing a static wall of mud to persist (Allouche, Frigaard & Sona Reference Allouche, Frigaard and Sona2000; Frigaard, Leimgruber & Scherzer Reference Frigaard, Leimgruber and Scherzer2003; Zare, Roustaei & Frigaard Reference Zare, Roustaei and Frigaard2017). These residual layers are called (wet) micro-annuli. As the cement hydrates, water is extracted from both mud layers and channels, potentially leaving behind a longitudinal pathway of porous dried mud. Aside from these, chemical shrinkage and debonding of cement from the walls produce (dry) micro-annuli. A combination of these defects is likely responsible for some of the larger leakage rates observed (Trudel & Frigaard Reference Trudel and Frigaard2023).
Fluid mechanics studies of primary cementing flows started in the 1960’s (McLean et al. Reference McLean, Manry and Whitaker1967). Although three-dimensional computational studies and lab-scale experimental studies are often performed (Malekmohammadi et al. Reference Malekmohammadi, Carrasco-Teja, Storey, Frigaard and Martinez2010; Renteria & Frigaard Reference Renteria and Frigaard2020; Sarmadi, Renteria & Frigaard Reference Sarmadi, Renteria and Frigaard2021; Zhang & Frigaard Reference Zhang and Frigaard2023), and give detailed information, the excessive length of cemented sections makes these approaches infeasible for understanding large-scale features of the process flows. Instead, the main workhorse for process design is a two-dimensional model that is based on averaging across the annular gap: effectively a Hele-Shaw cell with a non-uniform gap thickness (Bittleston, Ferguson & Frigaard Reference Bittleston, Ferguson and Frigaard2002). The Hele-Shaw, or two-dimensional gap-averaged (2DGA), approach was initially developed for laminar flows, assuming a uniform concentration of fluids across the annular gap. This style of model is useful in including both buoyancy and rheological effects, e.g. mud channels can be predicted. An iterative solution based on the augmented Lagrangian method is able to compute these nonlinear effects (Pelipenko & Frigaard Reference Pelipenko and Frigaard2004b ). Solutions show relevant process features, such as steady travelling wave displacement fronts that are also to some extent predictable theoretically (Pelipenko & Frigaard Reference Pelipenko and Frigaard2004a , Reference Pelipenko and Frigaardc ). This method has been extended to transitional, turbulent and mixed flow regimes, which happen commonly (Maleki & Frigaard Reference Maleki and Frigaard2017, Reference Maleki and Frigaard2018). These models have also been analysed for horizontal wells, where buoyancy effects are different (Carrasco-Teja et al. Reference Carrasco-Teja, Frigaard, Seymour and Storey2008), and have been extended to scenarios where the inner cylinder (casing) is moved (Carrasco-Teja & Frigaard Reference Carrasco-Teja and Frigaard2009, Reference Carrasco-Teja and Frigaard2010).
Models such as the 2DGA model have been used productively for many years, but have drawbacks. Inherent in the Hele-Shaw approach for laminar flows are assumptions that the streamlines are aligned in the flow direction and that inertial terms are small. Thus, multi-dimensional flow features (e.g. extreme rugosity of the outer wall, local effects near centralisers), are missed and there are many cementing flows where inertia becomes important in laminar regimes. This paper does not address these features.
The second and important limitation is that the fluid concentrations are not uniform across the gap thickness. Defining a Péclet number (
$Pe$
) based on the molecular diffusivity, a flow width
$\hat {D}$
and length
$\hat {L}$
, we find that
$Pe\gg {\hat {L}/\hat {D}}$
for most cementing flows. Thus, although fluids are miscible, fluids do not mix on the time scale of the displacement. Certainly we are far from the classical Taylor dispersion regime where
$1\ll Pe\ll {\hat {L}/\hat {D}}$
(Taylor Reference Taylor1953; Aris Reference Aris1956). The front between the two Newtonian fluids displacing in a channel remains very sharp, with the highest velocity at the centre of the channel, where the velocity is 1.5 times the mean velocity (Chen & Meiburg Reference Chen and Meiburg1996; Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Yang & Yortsos Reference Yang and Yortsos1997). This results in advective dispersion coming from the underlying Poiseuille velocity profile.
Yang & Yortsos (Reference Yang and Yortsos1997) developed asymptotic solutions for miscible Newtonian flows in a narrow plane channel to consider the impact of advective dispersion on the fluid flows. This approach was later extended by Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999) to include also buoyancy effects. In this approach, an assumed distribution of displaced and displacing fluids occurs across the channel. Asymptotic methods are then used to reduce the flow equations to a one-dimensional quasilinear conservation law that exhibits pure dispersion, kinematic shocks and other features that can also be found in experiments (Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999; Homsy, Scoffoni & Lajeunesse Reference Homsy, Scoffoni and Lajeunesse2001).
To bring this approach into the domain of cementing flows, Zhang & Frigaard (Reference Zhang and Frigaard2022) introduced a new approach, the so-called dispersive two-dimensional gap-averaged (D2DGA) model, to capture advective dispersion in a narrow eccentric annulus during the displacement of two Newtonian fluids. The novelty of this method comes from a two-layer flow assumption for the fluid concentration profile across the gap direction (outlined below). The comparison between results obtained from a three-dimensional simulation and numerical experiments with the D2DGA model showed a significant improvement compared with 2DGA simulations. Not only did this highlight the model’s accuracy, but also the necessity of including dispersion in modelling the displacement process. Favourable comparisons with experiments were made in Zhang & Frigaard (Reference Zhang and Frigaard2023).
While the D2DGA model of Zhang & Frigaard (Reference Zhang and Frigaard2022) brings a significant improvement over the classical 2DGA model, while retaining a two-dimensional formulation and consequent numerical efficiency, it neglects the influence of more complex rheological properties. This is needed for the purposes of simulating primary cementing flows more fully. Aside from model extension, there are underlying questions to answer and intuition that is missing once we consider dispersion. On the one hand, there is the intuitive notion that shear thinning and yield stress results in a blunter velocity profile, which may reduce dispersion. On the other hand, yield stress fluids can become stuck to the walls of a duct, which means that the mobile displacing fluid will move faster, increasing dispersion.
To build intuition and extend the D2DGA approach, this paper extends the D2DGA approach to pairs of Herschel–Bulkley fluids, which encompass the typical range of cementing fluids. This extension is a formidable modelling task, as although it is possible to write out the closure equations, their solution is not obvious. Even having developed the closure expressions, there are questions as to whether the two-dimensional streamfunction equation can be solved and how. Hence, here we only address a restricted planar displacement flow. Lastly comes the question of external verification, e.g. as was performed in Zhang & Frigaard (Reference Zhang and Frigaard2022, Reference Zhang and Frigaard2023), which is postponed to a future study. Therefore, what is presented is the development of a modelling framework.
An outline of the paper is as follows. The underlying principle of the D2DGA model is first reviewed in § 2.1. The impacts of the two-layer flow on the mobility and other closure functionals for Newtonian and non-Newtonian fluids are explained in §§ 2.2 and 2.3, respectively. The behaviour of a planar displacement front scenario is analysed in § 3 for wide ranges of rheological properties. In § 3.3 we classify the displacement process into (i) stable, (ii) partial penetration and (iii) unstable regimes using Muskat’s approach at various buoyancy numbers (
$b$
) and rheological parameters (
$n_k$
,
$\kappa _k$
). The paper ends with a brief summary in § 4.
2. Methodology
The 2DGA approach highlights the two-dimensional features of the flow by averaging the velocity and fluid concentration profiles in the radial direction. The derivation of the model is explained in detail in Bittleston et al. (Reference Bittleston, Ferguson and Frigaard2002), Maleki & Frigaard (Reference Maleki and Frigaard2017) and involves the use of standard scaling arguments to reduce the Navier–Stokes equations to the following dimensionless two-dimensional shear flow approximation:
\begin{align} 0&=-\frac {\partial p }{\partial y}, \end{align}
\begin{align} \displaystyle 0& = -\frac {1}{r_a}\frac {\partial p }{\partial \phi }+\frac {\partial \unicode{x03C4}_{\phi y}} {\partial y}+\frac {\rho \sin \beta \sin \pi \phi } {{\textit{Fr}}^{*^2}} , \end{align}
\begin{align} \displaystyle 0 &= -\frac {\partial p }{\partial \xi }+\frac {\partial \unicode{x03C4}_{\xi y}}{\partial y}-\frac {\rho \cos \beta } {{\textit{Fr}}^{*^2}}. \end{align}
Here
$\phi$
denotes the azimuthal coordinate and
$\xi$
measures distance along the axis of the inner pipe (casing);
$p$
is the pressure,
$\unicode{x03C4}_{\phi y}$
and
$\unicode{x03C4}_{\xi y}$
are the leading-order shear stresses,
$\rho$
is the dimensionless density and
$Fr^*$
is the Froude number defined by
\begin{equation} Fr^*=\sqrt {\frac {\hat {\unicode{x03C4}}_0} {\hat {\rho }_1 \hat {g} \hat {d}^*}}, \end{equation}
where
$\hat {\rho }_1$
is the density of fluid 1,
$\hat {g}$
is the gravitational acceleration and
$\hat {\unicode{x03C4}}_0$
represents a suitable stress scale for the flows. Throughout this paper we denote dimensional variables with the `
$\hat {\boldsymbol{\cdot }}$
’ accent and dimensionless variables without.
The inner and outer walls of the annulus have mean radii
$\hat {r}_i$
and
$\hat {r}_o$
, both of which might vary slowly with
$\xi$
. Axial and azimuthal distances have been scaled with
$\pi \hat {r}_a^*$
, where
$\hat {r}_a^*$
denotes a
$\xi$
-averaged mean radius. The radial distances are measured from a mean radial position and are scaled instead with a
$\xi$
-averaged half-annular gap width,
$\hat {d}^*$
, to define the scaled radial coordinate
$y$
; see figure 2 for a sketch of the geometry. In a typical cementing set-up a representative flow rate
$\hat {Q}_0$
can be defined for the process and used to define a mean velocity,
$\hat {Q}_0 = 4\pi \hat {r}_a^* \hat {d}^* \hat {w}_0$
. This, in turn, can be used with the rheological properties of the displaced fluid 1 to define
$\hat {\unicode{x03C4}}_0$
, e.g.
$\hat {\unicode{x03C4}}_0 = (\hat {\mu }_1\hat {\mu }_2)^{1/2} \hat {w}_0/\hat {d}^*$
for two Newtonian fluids. Equations (2.2) and (2.3) represent a two-dimensional shear flow in the plane of the annulus
$(\phi ,\xi )$
. After averaging across the annular gap and some simplifying assumptions, this model leads to a Hele-Shaw-type flow in a two-dimensional region of varying annular gap width. In general,
$y \in [-H,H]$
, where
$H=H(\phi ,\xi )$
is an
$O(1)$
function.
(a) Well geometry, (b) schematic of the eccentric annular space between the outer surface of the casing and the wellbore wall, (c) mapping the annulus to a Hele-Shaw cell.
A key assumption of the 2DGA model is that there is a single fluid present across the annular gap, i.e.
$\rho$
in (2.2) and (2.3) is a concentration-dependent (mixture) density and the constitutive laws also rely on concentration-dependent rheological parameters. This single fluid assumption has validity for turbulent shear flows (Maleki & Frigaard Reference Maleki and Frigaard2017), but in the laminar regime there will be a displacement front within the annular gap that evolves with the flow. Indeed such flows both exhibit dispersion and are vulnerable to various fingering and interfacial instabilities (Tabeling, Zocchi & Libchaber Reference Tabeling, Zocchi and Libchaber1987; Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996). In the primary cementing context, the fluids involved are viscous and non-Newtonian. The flows are characterised as either high-Péclet-number miscible or large-Capillary-number immiscible displacement flows. The latter are not considered here. Where displacement fronts evolve significantly in the flow directions, a modified version of the 2DGA model is needed to investigate dispersive effects. This modification has been made by Zhang & Frigaard (Reference Zhang and Frigaard2022) for two Newtonian fluids and extending the modelling approach for two Herschel–Bulkley fluids is the focus here.
2.1. Underlying principle of the dispersive two-layer model (D2DGA)
The mean annular gap,
$2\hat {d}=(\hat {r}_o-\hat {r}_i)$
, may vary slowly with
$\xi$
but remains of size
$\hat {d}^*$
. The narrow gap assumption is that this is small with respect to the mean circumference:
$2\pi \hat {r}^*_a$
, i.e.
$\hat {d}^*/(\pi \hat {r}^*_a)=\delta /\pi \ll 1$
. Generally speaking, the following three factors contribute to the 2DGA model limitations. (i) Elimination of inertial terms in (2.2) and (2.3), which formally needs
${Re}\delta /\pi \ll 1$
. (ii) Simplification of the flow to a shear flow, where terms of
$O(\delta /\pi )$
are neglected. (iii) Simplification of the mass transport equation to include only advective mass transport by the gap-averaged velocity. Both (i) and (ii) lead to approximations that are well understood (Bittleston et al. Reference Bittleston, Ferguson and Frigaard2002; Maleki & Frigaard Reference Maleki and Frigaard2017), but (iii) is less well explored. In the following, we first explain the underlying principle of the D2DGA model (Zhang & Frigaard Reference Zhang and Frigaard2022). Then, we extended this model to two viscoplastic fluids.
With the same lubrication theory scaling assumptions as before, the main difference to the 2DGA model is in the assumed mass transport. Considering a high-Péclet-number (
$Pe$
) miscible laminar flow, the thickness of diffusive layers scales with
$\sqrt {1/ Pe}$
, meaning that distinct fluid layers may persist across the annulus gap. To capture the impact of advective dispersion on the displacement, it is assumed that fluid 2 can disperse symmetrically along the centre of the annular gap, between interfaces at
$y=\pm y_i$
, and a thin film of fluid 1 remains near the walls, in
$[-H, y_i)$
and
$(y_i,H]$
, as shown in figure 3. Thus, in the D2DGA model a jump of fluid properties at the interface appears, which is neglected in the 2DGA model via the assumption of transverse mixing. The dimensionless governing equations for fluid 1 (displaced fluid) and fluid 2 (displacing fluid) are
\begin{equation} \mbox{Fluid 1} \hspace {0.1 cm} \begin{cases} \displaystyle 0=-\frac {1}{r_a}\frac {\partial p }{\partial \phi }+\frac {\partial \unicode{x03C4}_{\phi y, 1}} {\partial y}+\frac {\rho _1 \sin \beta \sin \pi \phi } {{\textit{Fr}}^{*^2}}, \\[10pt] \displaystyle 0=-\frac {\partial p }{\partial \xi }+\frac {\partial \unicode{x03C4}_{\xi y, 1}} {\partial y}-\frac {\rho _1 \cos \beta } {{\textit{Fr}}^{*^2}}, \quad y\in ( y_i, H], \end{cases} \end{equation}
\begin{align} \mbox{Fluid 2} \hspace {0.1 cm} \begin{cases} \displaystyle 0=-\frac {1}{r_a}\frac {\partial p }{\partial \phi }+\frac {\partial \unicode{x03C4}_{\phi y, 2}} {\partial y}+\frac {\rho _2 \sin \beta \sin \pi \phi } {{\textit{Fr}}^{*^2}}, \\ \displaystyle 0=-\frac {\partial p }{\partial \xi }+\frac {\partial \unicode{x03C4}_{\xi y, 2}} {\partial y}-\frac {\rho _2 \cos \beta } {{\textit{Fr}}^{*^2}}, \quad y\in [0, y_i). \end{cases} \\[6pt] \nonumber \end{align}
Schematic of 2DGA and D2DGA models.
Note that in (2.5) and (2.6) the equations are valid only within each fluid layer, as indicated, but the densities and rheological properties are also constant in each layer. In addition, the radial momentum (2.1) is satisfied across each layer. We now introduce the density difference between the two fluids,
$\Delta \hat {\rho } =\hat {\rho }_1 -\hat {\rho }_2$
, and the mean density,
\begin{align} \hat {\rho }=(1-\bar {c})\hat {\rho }_1 +\bar {c}\hat {\rho }_2 = \hat {\rho }_1 - \bar {c} \Delta \hat {\rho }, \end{align}
as well as their dimensionless versions, obtained by dividing with
$\hat {\rho }_1$
. We define the concentration
$\bar {c}= y_i / H$
such that it corresponds directly to the fraction of displacing fluid; see figure 3.
The leading-order velocity is
$(v,w)$
in
$(\phi ,\xi )$
directions, and the leading-order strain rate comes from the
$y$
derivatives of
$(v,w)$
. The derivatives of
$(v,w)$
in
$(\phi ,\xi )$
directions come from variations in the annular gap geometry, which are small due to scaling arguments. For a generalised Newtonian fluid
$k$
,
\begin{equation} (\unicode{x03C4}_{\phi y, k},\unicode{x03C4}_{\xi y, k}) = \eta _k(\dot {\gamma }) \left ( \frac {\partial v}{\partial y}, \frac {\partial w} {\partial y} \right ),\,\,\,\,\,\,\dot {\gamma } = \left [ \left ( \frac {\partial v}{\partial y}\right )^2 + \left ( \frac {\partial w} {\partial y} \right )^2 \right ]^{1/2}, \end{equation}
where
$\tau _k = \eta _k(\dot {\gamma })\dot {\gamma }$
is the constitutive law and
$\tau _k = [\unicode{x03C4}_{\phi y, k}^2,\unicode{x03C4}_{\xi y, k}^2 ]^{1/2}$
. Note that for monotone flow curves, we may equally write the effective viscosity as
$\eta _k =\eta _k(\tau _k)$
, or below as
$\eta _k =\eta _k(y)$
, implicitly recognising the linear stress variation with
$y$
.
By integrating (2.5) and (2.6) with respect to
$y$
we arrive at
\begin{align}& (\unicode{x03C4}_{\phi y, 2},\unicode{x03C4}_{\xi y, 2}) = \boldsymbol{\unicode{x03C4}_2}=[-\boldsymbol{G}+(1-\bar {c})\boldsymbol{G_b}]y,\quad \hspace {1 cm} y\in [0, y_i], \end{align}
\begin{align}& (\unicode{x03C4}_{\phi y, 1},\unicode{x03C4}_{\xi y, 1}) = \boldsymbol{\unicode{x03C4}_1} =\boldsymbol{\unicode{x03C4}_i}-[\boldsymbol{G}+\bar {c}\boldsymbol{G_b}](y-y_i), \hspace {1 cm} y\in (y_i, H]. \\[9pt] \nonumber \end{align}
Here
$\boldsymbol{\unicode{x03C4}_i}$
is the interfacial shear stress,
$\boldsymbol{\unicode{x03C4}_i}=[-\boldsymbol{G}+(1-\bar {c})\boldsymbol{G_b}]y_i$
. The expressions for
$\boldsymbol{G}$
and
$\boldsymbol{G_b}$
are defined as
\begin{align}& \boldsymbol{G}=-\left(\frac {1}{r_a}\frac {\partial p }{\partial \phi }-\frac {\rho \sin \beta \sin \pi \phi } {{\textit{Fr}}^{*^2}},\frac {\partial p }{\partial \xi }+\frac {\rho \cos \beta } {{\textit{Fr}}^{*^2}}\right), \end{align}
\begin{align}&\qquad \boldsymbol{G_b}= - b (f_\xi , -f_\phi ): \,\,\,\, b = -\frac {\Delta \rho }{r_a {\textit{Fr}}^{*^2}}, \end{align}
where
$\boldsymbol{f}=(f_\phi ,f_\xi )= (r_a \cos \beta ,r_a \sin \beta \sin \pi \phi )$
. The vector
$\boldsymbol{G}$
we shall call the modified pressure gradient. Observe that in the case of the 2DGA model,
$\boldsymbol{\unicode{x03C4}} = -\boldsymbol{G}y$
(see (2.2) and (2.3)), i.e.
$\boldsymbol{G}$
is the frictional pressure gradient felt by the volume-averaged mixture of the two fluids. However, for the D2DGA model, there is also the vector
$\boldsymbol{G_b}$
that explicitly represents the buoyancy force. Note that at
$y=H$
the buoyancy term (
$\boldsymbol{G_b}$
) plays no part and
$\boldsymbol{\unicode{x03C4}_1}(H)=-\boldsymbol{G}H$
, i.e.
$|\boldsymbol{G}|H$
gives the wall shear stress, just as in the 2DGA model.
We have separated
$\boldsymbol{G_b}$
into its direction and amplitude. Regarding the direction, note that both
$r_a$
and
$\beta$
are assumed to vary only slowly with
$\xi$
. The consequence of this is that primarily it is gradients in concentration that produce buoyancy effects and not changing orientation of the well. The amplitude
$b$
has been signed so that typically we have
$b\gt 0$
in a primary cementing flow, where denser fluids typically displace lighter fluids in the upwards direction. Observe that
$b$
represents a ratio of the buoyancy stress (
$-\Delta \hat {\rho } \hat {g} \hat {d}^*$
) to the viscous stress used to scale the flow (
$\hat {\unicode{x03C4}}_0$
).
Equation (2.10) indicates that the jump in fluid properties at the interface between the two fluids leads to a deviation in shear stress direction within the displaced fluid, compared with that of the displacing fluid. Since the shear stresses determine the strain rate, via constitutive laws for each fluid, this directional jump will also be felt in the fluid velocity components. To elaborate more, the buoyancy vector can independently change the flow direction, and the mean flow is not along the direction of the modified pressure gradient (
$\boldsymbol{G}$
). In figure 4(a) as an example, we have plotted the stress components that result from (2.9) and (2.10) for a particular parameter set. These equations have used symmetry at
$y=0$
and continuity of shear stress at
$y=y_i$
. The distance across the channel is indicated by the varying colour of the lines. We see the two linear segments, representing the two fluid layers, for the case of a non-zero
$\boldsymbol{G_b}$
. The single linear line comes from using only
$\boldsymbol{G}$
, (with
$\boldsymbol{G_b}=0$
), which can also be thought of as representing the stress variation in the 2DGA model. At the wall the two realisations coincide, since the wall shear stress is determined only by
$\boldsymbol{G}H$
.
(a) Variation of shear stresses at
$\boldsymbol{G}=[6\hspace {0.1 cm}1]$
,
$\boldsymbol{n}=[1 \hspace {0.1cm}1]$
,
$\boldsymbol{\unicode{x03C4}_Y}=[0 \hspace {0.1cm}0]$
,
$\boldsymbol{\kappa }=[1 \hspace {0.1cm}1], \bar {c}=0.5$
, and (b) shear strain at
$\boldsymbol{G}=[6\hspace {0.1 cm}1]$
,
$\boldsymbol{n}=[0.5 \hspace {0.1cm}0.5]$
,
$\boldsymbol{\unicode{x03C4}_Y}=[1 \hspace {0.1cm}1]$
,
$\boldsymbol{\kappa }=[1 \hspace {0.1cm}1], \bar {c}=0.5$
for various
$0\leqslant y\leqslant 1$
, the colour bar illustrates values associated to
$y$
. In the straight line
$\boldsymbol{G_b}=[0 \hspace {0.1cm}0]$
and in the deviated line
$\boldsymbol{G_b}=[1 \hspace {0.1cm}1]$
. The change in the direction of flow comes from the sudden jump in fluid properties at the interface (
$\boldsymbol{G_b} y_i/y$
).
Using the constitutive laws of a Herschel–Bulkley fluid for each fluid, we may also compute the two components of the velocity gradient as we traverse the channel. These are shown in figure 4(b). Again two examples are shown: with and without non-zero
$\boldsymbol{G_b}$
. Here the fluid properties selected include a yield stress in both fluids. Due to the symmetry condition at
$y=0$
, the yield stress is not exceeded close to the centre of the channel and we see that a finite interval of
$y$
values map to zero velocity gradient, i.e. a central plug region, for both cases. Again the case with non-zero buoyancy shows a jump in velocity gradients at the interface, whereas the isodense case keeps a constant ratio between the two velocity gradients.
From the velocity gradient, using the no slip boundary condition at
$y=H$
, and continuity of velocity at
$y=y_i$
, the velocity profile may be calculated directly (see Appendix A). Instead, the gap-averaged velocity can be calculated directly, using integration by parts, as in Zhang & Frigaard (Reference Zhang and Frigaard2022):
\begin{equation} H\boldsymbol{\bar {u}}=(H\bar {v}, H\bar {w})=\mathcal{I}_1\boldsymbol{G}-\frac {\mathcal{I}_2}{H}\boldsymbol{G_b}. \end{equation}
The functions
$\mathcal{I}_1$
and
$\mathcal{I}_2$
are
\begin{align}&\qquad\quad \mathcal{I}_1=\int _{0}^{y_i} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y}+ \int _{y_i}^{H} \frac {\tilde {y}^2}{\eta _1(\tilde {y})} {\rm d}\tilde {y} = H^3 I_1(\bar {c};\boldsymbol{\cdot }) , \end{align}
\begin{align}& \mathcal{I}_2=(H-y_i)\int _{0}^{y_i} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y}+y_i\int _{y_i}^{H} \frac {\tilde {y}(H-\tilde {y})}{\eta _1(\tilde {y})} {\rm d}\tilde {y} = H^4 I_2(\bar {c};\boldsymbol{\cdot }) . \end{align}
Note that, unlike the Newtonian case in Zhang & Frigaard (Reference Zhang and Frigaard2022), the effective viscosity of the fluids varies with the shear stress, which means that (2.14) and (2.15) also depend explicitly on
$\boldsymbol{G}$
and
$\boldsymbol{G_b}$
. The expressions for
$ I_1(\bar {c};\boldsymbol{\cdot })$
and
$ I_2(\bar {c};\boldsymbol{\cdot })$
are given later in Appendix A. Unlike Zhang & Frigaard (Reference Zhang and Frigaard2022), here we use
$ I_1$
and
$ I_2$
in place of
$\mathcal{I}_1$
and
$\mathcal{I}_2$
as they separate out the dependency on the local gap thickness from that on
$c$
and the dimensionless rheological groups, represented above by the `
$\boldsymbol{\cdot}$
’. Observe that physically both
$I_1$
and
$I_2$
represent mobilities. We refer to
$I_1$
as the mean mobility and
$I_2$
as the buoyant mobility.
The streamfunction for the incompressible flow is defined as
${\nabla }_a\varPsi =(2H\bar {w},-2H\bar {v})$
. By rearranging (2.13) based on the streamfunction and cross-differentiating, we eliminate the pressure and obtain the partial differential equation
\begin{equation} {\nabla }_a \boldsymbol{\cdot }[\boldsymbol{S}+\boldsymbol{b}]=0, \end{equation}
which is analogous to the 2DGA model. Here
$\boldsymbol{S}$
contains pressure gradient terms and is defined as
\begin{equation} \boldsymbol{S} = \frac {r_a {\nabla }_a\varPsi }{2\mathcal{I}_1} = \frac {r_a {\nabla }_a\varPsi }{2H^3 I_1}. \end{equation}
Assuming that the flow is non-zero it can be shown that there is a unique mapping between
$\boldsymbol{G}$
and
${\nabla }_a\varPsi$
(see Appendix A), so that
$\mathcal{I}_1$
and
$\mathcal{I}_2$
may instead be considered to depend implicitly on
${\nabla }_a\varPsi$
.
The second term, in the direction
$\boldsymbol{f}$
, can be written as
\begin{equation} \boldsymbol{b} = \left (\frac {\rho }{{\textit{Fr}}^{*^2}}+\frac {b}{H}\frac {\mathcal{I}_2}{\mathcal{I}_1}\right ) \boldsymbol{f} = \left (\frac {1}{{\textit{Fr}}^{*^2}}+b\left [\bar {c} + \frac {I_2}{I_1} \right ]\right ) \boldsymbol{f}. \end{equation}
Observe that with our assumption on the slow variation with
$\xi$
, the first term can be neglected, i.e. since
${\nabla }_a \boldsymbol{\cdot }\boldsymbol{f} = 0$
. The second term
$b \bar {c} \boldsymbol{f}$
gives the buoyancy gradients of the mean density and is also present in the 2DGA model. The last part comes directly from the layered flow assumption, i.e. advective dispersion directly gives rise to a buoyancy source term (
$ b \boldsymbol{f} \boldsymbol{\cdot }{\nabla }_a [I_2/I_1]$
), that is completely absent from the 2DGA model.
2.1.1. Evolution of the concentration
To advance the fluid concentration, the gap-averaged velocity in the transport equation is replaced by a flux function
\begin{equation} \frac {\partial [H r_a \bar {c}]}{\partial t }+ {\nabla }_a \boldsymbol{\cdot }\boldsymbol{q}=0, \end{equation}
where
$\boldsymbol{q}$
is the volumetric flux passing through the displacing fluid layer:
\begin{equation} \boldsymbol{q}=(q_\phi ,q_\xi )=r_a\left(\int _{0}^{y_i} v {\rm d}y,\int _{0}^{y_i} w {\rm d}y\right)\!. \end{equation}
The volumetric flux function can be evaluated by expressing
$(v,w)$
in terms of integrals of their
$y$
derivatives, then using (2.9) and (2.10), and finally, (2.13) to eliminate the modified pressure gradient. The expression is
\begin{equation} (q_\phi ,q_\xi )=\frac {r_a}{2}\left(-\frac {\partial \varPsi }{\partial \xi }, \frac {1}{r_a}\frac {\partial \varPsi }{\partial \phi }\right)q_0(\bar {c}) - b H^3 I_3(\bar {c}) (f_\xi ,-f_\phi ). \end{equation}
The isotropic flux
$q_0(\bar {c})$
is given by
\begin{eqnarray} q_0(\bar {c};\boldsymbol{\cdot }) & = & \frac {1}{H^3 I_1} \displaystyle { \left ( \int _{0}^{y_i} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y} + y_i\int _{y_i}^{H} \frac {\tilde {y}}{\eta _1(\tilde {y})} {\rm d}\tilde {y} \right ) }. \end{eqnarray}
By construction,
$q_0(0;\boldsymbol{\cdot }) = 0$
and
$q_0(1;\boldsymbol{\cdot }) =1$
. Note that the
$H^3$
power normalises the two integrals. In the case of the 2DGA model,
$q_0(\bar {c};\boldsymbol{\cdot }) = \bar {c}$
. The buoyancy component of the volumetric flux is governed by
$I_3(\bar {c};\boldsymbol{\cdot })$
:
\begin{align} I_3(\bar {c};\boldsymbol{\cdot }) &= \frac {1}{H^7I_1} \left [ \left (\int _{0}^{y_i} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y}+y_i\ \int _{y_i}^{H} \frac {\tilde {y}}{\eta _1(\tilde {y})} {\rm d}\tilde {y} \right ) \mathcal{I}_2 \right . \nonumber \\ & \quad \left . - \left ( (H-y_i) \int _{0}^{y_i} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y}+y_i^2\int _{y_i}^{H} \frac {(H-\tilde {y})}{\eta _1(\tilde {y})} {\rm d}\tilde {y} \right ) \mathcal{I}_1 \right ] \nonumber \\ & = \frac {1}{H^7I_1}\left [ \left ( y_i\int _{0}^{y_i} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y} \int _{y_i}^{H} \frac {(\tilde {y}-y_i)}{\eta _1(\tilde {y})}{\rm d}\tilde {y} \right ) + \right . \nonumber \\ & \quad \left . + y_i^2 \left (\left [ \int _{y_i}^{H} \frac {\tilde {y}}{\eta _1(\tilde {y})} {\rm d}\tilde {y} \right ]^2 - \left [ \int _{y_i}^{H} \frac {1}{\eta _1(\tilde {y})} {\rm d}\tilde {y} \right ]\left [ \int _{y_i}^{H} \frac {\tilde {y}^2}{\eta _1(\tilde {y})} {\rm d}\tilde {y} \right ] \right ) \right ] .\end{align}
Again, the
$H^7$
power normalises all the integral terms, which now depend only on
$c$
and rheological parameters. We refer to
$I_3$
as the buoyancy flux distribution.
2.2. Newtonian fluids
Zhang & Frigaard (Reference Zhang and Frigaard2022) evaluated
$\mathcal{I}_1$
and
$\mathcal{I}_2$
analytically for two Newtonian fluids. Our formulae are slightly different, as our scaling is chosen with a view to the non-Newtonian fluid pairs, considered later. We use the geometric mean viscosity
$(\hat {\mu }_1\hat {\mu }_2)^{0.5}$
to define the stress scale
$\hat {\unicode{x03C4}}_0$
. The viscosity ratio is defined as
$m = \hat {\mu }_1/\hat {\mu }_2$
. Hence, we find that
$\eta _1 = m^{0.5}$
,
$\eta _2 = m^{-0.5}$
and the following expressions result:
\begin{align} \mathcal{I}_1 &= \left(\frac {y_i^3}{3\eta _2}+\frac {H^3-y_i^3}{3\eta _1}\right) = H^3 \frac {m^{0.5}\bar {c}^3 + m^{-0.5}\left(1-\bar {c}^3\right)}{3} = H^3 I_1(\bar {c};m), \end{align}
\begin{align} \mathcal{I}_2&=\left[\frac {(H-y_i)y_i^3}{3\eta _2}+\frac {y_i(H-y_i)^2(H+2y_i)}{6\eta _1}\right] \nonumber \\ &= H^4 \frac {2m^{0.5}\bar {c}^3(1-\bar {c})+m^{-0.5}\bar {c}(1-\bar {c})^2(1+2\bar {c})}{6} = H^4 I_2(\bar {c};m). \end{align}
The flux function in (2.21) can be simplified to
\begin{equation} q_0(\bar {c};m) = \bar {c}\left[\frac {m\bar {c}^2+1.5(1-\bar {c}^2)}{m\bar {c}^3+1-\bar {c}^3}\right] \end{equation}
and
\begin{equation} I_3(\bar {c};m)=\frac {\bar {c}^2(1-\bar {c})^3\left[4m\bar {c}+3(1-\bar {c}^2)\right]}{12[m\bar {c}^3+1-\bar {c}^3]}. \end{equation}
Observe that, for
$m=1$
and
$\Delta \rho =0$
, the flux function is equal to
\begin{equation} (q_\phi ,q_\xi )=\frac {r_a}{2}\left(-\frac {\partial \varPsi }{\partial \xi }, \frac {1}{r_a}\frac {\partial \varPsi }{\partial \phi }\right)\bar {c}(1.5-0.5\bar {c}^2), \end{equation}
which is
$(1.5-0.5\bar {c}^2)$
times the flux function coming from the 2DGA model, where
$q_0 = \bar {c}$
. This amplified term is associated with the underlying Poiseuille flow, which disperses the flow. This is similar to the flux term previously calculated analytically using a transverse flow equilibrium assumption to capture the dispersion of two miscible Newtonian fluids in a long channel (Yang & Yortsos Reference Yang and Yortsos1997). Note however that here the gap thickness (
$H$
) in an eccentric annulus depends on both azimuthal and axial directions, and the local mean velocity may also vary.
Variation of
${I}_1(\bar {c})$
and
$q_0(\bar {c})$
by
$\bar {c}$
for Newtonian fluids at
$m=5,\,1,\,0.2$
,
${\bar {w}}=1$
,
${b}=0 \text{ (isodensity)}$
.
Comparing the mean mobility
$I_1$
expression from the 2DGA model (
$I_1= [m^{0.5}\bar {c} + m^{-0.5}(1-\bar {c})] / 3$
) with (2.24), we can interpret the D2DGA model as a form of interpolation of the viscosity at intermediate concentrations; see figure 5(a). In the absence of the buoyancy force,
$q_0(\bar {c})$
determines the frontal behaviour of the displacement; see figure 5(b). Since
$d^2q_0 /{\rm d}\bar {c}^2\lt 0$
for the D2DGA model, the front with lower concentration travels at the higher speed, which disperses the interface. The straight line in figure 5(b) corresponds to the uniform front speed of the 2DGA model. This clearly highlights the limitation of the 2DGA model in capturing the dispersive fronts. Evidently, prediction of non-dispersive frontal behaviours using the 2DGA model can result in overestimating the volumetric efficiency, as noted in Zhang & Frigaard (Reference Zhang and Frigaard2022).
The effect of the viscosity ratio on the variation of mean mobility
$I_1(\bar {c})$
and isotropic flux
$q_0(\bar {c})$
are illustrated in figures 5(a) and 5(b) for the D2DGA model. The constant value of mean mobility for
$m=1$
is the same as that from the 2DGA model. On the other hand, whenever
$m\neq 1$
, the new expression defined for the mean mobility of fluids (
$I_1(\bar {c})$
) leads to a variation with concentration, between the 2DGA values, which are recovered at
$\bar {c} = 0$
and
$\bar {c} = H$
. When the viscosity of the displacing fluid is higher than the displaced fluid (
$m\lt 1$
), the front becomes more stable and the mean mobility of fluids decreases with
$\bar {c}$
. For invasion of a low viscous fluid into a high viscous fluid (
$m\gt 1$
), the front is less stable and the mean mobility of fluids increases with
$\bar {c}$
. We can also observe this in the behaviour of the isotropic flux
$q_0(\bar {c})$
. The wavespeeds (
$q_0^\prime (\bar {c})$
) at low
$\bar {c}$
increase with
$m$
, implying increased dispersion of the finger of fluid 2 into fluid 1. The buoyant mobility component
$I_2$
always vanishes at both
$\bar {c}=0$
and
$1$
, since buoyancy arises only from having two fluids present, but here this is irrelevant as
$b=0$
.
2.3. Non-Newtonian fluids
Cementing fluids are typically modelled as Herschel–Bulkley fluids, in which the consistency, power-law index and yield stress of fluid
$k$
are denoted by
$\hat {\kappa }_k$
,
$n_k$
and
$\hat {\unicode{x03C4}}_{Y,k}$
, respectively:
$k=1$
for the displaced fluid and
$k=2$
for the displacing fluid. As an effective viscosity scale we use
\begin{equation} \hat {\mu }_{e}=\left[\hat {\kappa }_{1}\hat {\dot {\gamma }}_0 ^{n_1-1}\hat {\kappa }_{2}\hat {\dot {\gamma }}_0 ^{n_2-1}\right]^{0.5}, \qquad \hat {\dot {\gamma }}_0 = \hat {w}_0/\hat {d}^* . \end{equation}
We then define the stress scale
$\hat {\unicode{x03C4}}_{0}$
as
\begin{equation} \hat {\unicode{x03C4}}_{0} = \hat {\mu }_{e} \hat {\dot {\gamma }}_0 + \max \{\hat {\unicode{x03C4}}_{Y,1},\hat {\unicode{x03C4}}_{Y,2} \}. \end{equation}
The scaled constitutive laws for fluid
$k$
are
\begin{equation} \begin{array}{lcl} \boldsymbol{\unicode{x03C4}_k} = \displaystyle {\left ( \kappa _k \left |\frac {{\rm d}\boldsymbol{u}}{{\rm d}y}\right |^{n_k-1} + \unicode{x03C4}_{Y,k} \left |\frac {{\rm d}\boldsymbol{u}}{{\rm d}y}\right |^{-1} \right ) \frac {{\rm d}\boldsymbol{u}}{{\rm d}y}}, &\,\,\,\,& \Leftrightarrow |\boldsymbol{\unicode{x03C4}_k}| \gt \unicode{x03C4}_{Y,k}, \\[1ex] \displaystyle {\frac {{\rm d}\boldsymbol{u}}{{\rm d}y}} = 0,& & \Leftrightarrow |\boldsymbol{\unicode{x03C4}_k}| \leqslant \unicode{x03C4}_{Y,k}. \end{array} \end{equation}
The scaled parameters can be written as
\begin{equation} \kappa _1 = \frac {\hat {\kappa }_{1}\hat {\dot {\gamma }}_0 ^{n_1}}{\hat {\unicode{x03C4}}_{0}} = \frac {m^{0.5}}{1+B}, \qquad \kappa _2 = \frac {\hat {\kappa }_{2}\hat {\dot {\gamma }}_0 ^{n_2}}{\hat {\unicode{x03C4}}_{0}} = \frac {m^{-0.5}}{1+B}, \qquad\unicode{x03C4}_{Y,k} = \frac {\hat {\unicode{x03C4}}_{Y,k}}{\hat {\unicode{x03C4}}_{0}} , \end{equation}
where
$m = \hat {\kappa }_{1}\hat {\dot {\gamma }}_0 ^{n_1-n_2}/\hat {\kappa }_{2}$
is the generalisation of the Newtonian viscosity ratio. The Bingham number
$B$
here is
\begin{align} B = \frac { \max \{\hat {\unicode{x03C4}}_{Y,1},\hat {\unicode{x03C4}}_{Y,2} \} }{\hat {\mu }_{e} \hat {\dot {\gamma }}_0} . \end{align}
Note that the parameters above simplify to the Newtonian definitions. The point of using
$\max \{\hat {\unicode{x03C4}}_{Y,1},\hat {\unicode{x03C4}}_{Y,2} \}$
is that the resulting scaled yield stresses are bounded within
$[0,1)$
, which is useful for numerical evaluation and tabulation.
Note that each
$I_k(\bar {c})$
and
$q_0(\bar {c})$
now depend on the constitutive laws of the fluids. In yielded regions,
$\eta _k(\tilde {y})$
is defined as
\begin{equation} \begin{array}{lcl} { \eta _k = \displaystyle {\left ( \kappa _k \left |\frac {{\rm d}\boldsymbol{u}}{{\rm d}y}\right |^{n_k-1} + \unicode{x03C4}_{Y,k} \left |\frac {{\rm d}\boldsymbol{u}}{{\rm d}y}\right |^{-1} \right ) }}. \end{array} \end{equation}
It follows that (2.13) becomes nonlinear, either considered as an equation in terms of
$\boldsymbol{\bar {u}}$
(equivalently
${\nabla }_a \varPsi$
) or
$\boldsymbol{G}$
. An explicit analytical solution is not possible in general, and the solution
$\boldsymbol{u}$
is needed in order to evaluate the different closure functionals (
$I_k(\bar {c})$
and
$q_0(\bar {c})$
). The solution can be found, either for fixed
$\boldsymbol{\bar {u}}$
or fixed
$\boldsymbol{G}$
, via the iterative solution of an equivalent convex minimisation problem. This is discussed in Appendix A.
As an example, figure 6 displays the variation of the mean mobility and isotropic flux for two shear-thinning fluids with equal power-law indexes and with no buoyancy. As the effective viscosities of the two fluids are identical there is no
$\bar {c}$
dependency in figure 6(a). As the power-law index of both fluids reduces (
$n\to 0.2$
), the effective viscosity decreases and leads to an increase in the mean mobility value. The variation of isotropic flux is illustrated in figure 6(b). The broken straight line corresponds to a planar interface between the two fluids. At the lowest power-law index of the fluids, the isotropic flux gets closer to the broken straight line. This coincides with the fluid velocity becoming progressively plug-like, as observed in the inset of figure 6(b).
Variation of mean mobility
$I_1(\bar {c})$
and isotropic flux and
$q_0(\bar {c})$
with
$\bar {c}$
based on D2DGA model for shear-thinning fluids at
$m=1, \bar {w}=1, b=0 \text{ (isodensity)}.$
As the power-law index of the fluids reduces, the effective viscosity
$(\hat {\mu }_e)$
decreases leading to a growth in the mean mobility value. The first and second components of
$\boldsymbol{n}=[n_1 \,n_2]$
represent the power-law index of the displaced and displacing fluids, respectively.
The impact of yield stress of both fluids on mean mobility and isotropic flux are shown in figure 7, again with no buoyancy. In figure 7(a,b) the yield stress of the displaced fluid is fixed at
${\unicode{x03C4}}_{Y,1}=0.94$
, and the yield stress of the displacing fluid is varied:
${\unicode{x03C4}}_{Y,2}=0.23,\, 0.33, \,0.71, \,0.94$
. Looking at the low range of
$\bar {c}$
in figure 7(a), the yield stress of the displacing fluid is not exceeded. The contribution of the displacing fluid (i.e.
$\int _{0}^{y_i} [\tilde {y}^2/\eta _2(\tilde {y})]\,{\rm d}\tilde {y}$
) to the mean mobility is therefore zero since
$\eta _2 \to \infty$
and
$\dot {\gamma }\to 0$
. Therefore, the mean mobility of the fluids is initially constant and comes only from the displaced fluid (i.e.
$\int _{y_i}^{H} [\tilde {y}^2/\eta _1(\tilde {y})]\,{\rm d}\tilde {y}$
), and specifically from a yielded layer close to the wall. This contribution is identical to that for when the channel is fully filled with the displaced fluid
$I_1(\bar {c})\approx 0.9$
. To illustrate these constant regions at small
$\bar {c}$
, we have separated out the contributions to the mean mobility (
$I_1$
) from displaced and displacing fluids and plotted these against
$\bar {c}$
in the inset of figure 7(a) for
$\boldsymbol{\tau _Y}=[0.94\,0.23], \boldsymbol{n}=[0.5\,0.5]$
.
The effect of the yield stress of the displacing fluid (a,b) and the displaced fluid (c,d) on the variation of mean mobility
$I_1(\bar {c})$
and isotropic flux
$q_0(\bar {c})$
with
$\bar {c}$
. The results are obtained using the D2DGA model for Herschel–Bulkley fluids and isodense flows (
$b=0$
) at
$m=1, B=18.3$
for
$n=[0.5 \hspace {0.1 cm}0.5]$
and
$B=17.66$
for
$n=[0.7 \hspace {0.1 cm}0.7$
]. The first and second components of
$\boldsymbol{n}=[n_1 \,n_2]$
and
$\boldsymbol{\tau }_Y=[\tau _{Y,1} \,\tau _{Y,2}]$
represent the power-law index and the yield stress of the displaced and displacing fluids, respectively.
As the concentration increases, the displacing fluid starts yielding, and the impact of the mobility of the displacing fluid on the mean mobility (
$I_1(\bar {c})$
) increases. We observe that
$I_1$
is approximately linear over a range of
$\bar {c}$
, with growth rate
$a$
, which reduces as the yield stress of the displacing fluid approaches that of the displaced fluid (
${\unicode{x03C4}}_{Y,2}\to 0.94$
). From the inset of figure 7(a), we see that the approximately linear growth of the mean mobility starts as the concentration exceeds
$\bar {c}\approx 0.52$
, and coincides with the emergence of a static wall layer in fluid 1 (within which the mobility of fluid 1 is zero).
Figure 7(b) displays the impact of the yield stress of the displacing fluid on the frontal behaviour, felt via the isotropic flux function
$q_0(\bar {c})$
. The rheological parameters are identical with figure 7(a) and the fluids have the same density (
$b=0$
). At small values of
$\tau _{Y,2}$
, e.g.
$0.23,\,0.33$
, due to the high yield stress of the displaced fluid, there is an interval of
$\bar {c}$
for which the displaced fluid is completely unyielded and, consequently, the isotropic flux increases linearly with
$\bar {c}$
(given by the plug velocity multiplying
$\bar {c}$
). As the displacing fluid begins to yield, we see a transfer to fluid 2. The minimum concentration for which
$q_0(\bar {c}) = 1$
is denoted
$\bar {c}_{\textit{min}}$
and can be directly calculated. For
$\bar {c} \gt \bar {c}_{\textit{min}}$
, fluid 1 remains static and
$q_0(\bar {c}) = 1$
.
In figure 7(c,d) the impact of the yield stress of the displaced fluid is investigated, where the yield stress of the displacing fluid is fixed at
${\unicode{x03C4}}_{Y,2}=0.94$
. Unlike Figure7(a), the mean mobility
$I_1$
is approximately
$\bar {c}$
independent for a wide range of concentrations. Fluids with lower
$\tau _{Y,1}$
have higher mean mobility values. Figure 7(d) shows the analogous variations in isotropic flux, which are close to linear. This comes from the high plasticity and central plug region of the displacing fluid.
2.3.1. Variations of mean mobility (
$I_1(\bar {c})$
) with static wall layers
We have still not explained the apparent linear growth of
$I_1$
for
$\bar {c}\geqslant \bar {c}_{\textit{min}}$
, although this appears to be associated with the regime of static wall layers. Assuming that
$\bar {c}\geqslant \bar {c}_{\textit{min}}$
, only fluid 2 is flowing. Since the flow is still computed for a fixed flow rate, but now with a single fluid, we may use the earlier results from the 2DGA model to find the relation between
$G_\xi$
and
$\bar {c}$
, i.e. the displacing fluid is injected into a channel with a half-width
$\bar {c}$
at a mean velocity of
$\bar {w}=1$
:
\begin{align} \int _{0}^{\bar {c}} w {\rm d}\tilde {y}=F_2(G_\xi ). \end{align}
Here
$F_2(G_\xi )$
is defined as
\begin{equation} F_2(G_\xi )= \displaystyle \frac {\left(\left(\dfrac {1}{n_2}+1\right)\bar {c}G_\xi +\unicode{x03C4}_{Y,2}\right)}{G^2_\xi \left(\dfrac {1}{n_2}+1\right)\left(\dfrac {1}{n_2}+2\right)}\dfrac {(\bar {c}G_\xi -\unicode{x03C4}_{Y,2})_{+}^{\left(\frac {1}{n_2}+1\right)}}{\kappa ^{(1/n_2)}}, \quad \bar {c}G_\xi \gt \unicode{x03C4}_{Y,2}. \end{equation}
When the displaced fluid is static, (2.14) can also be simplified to (2.37) using integration by parts, such that the flow rate constraint is satisfied:
\begin{equation} {I}_1(\bar {c};.)=\int _{0}^{\bar {c}} \frac {\tilde {y}^2}{\eta _2(\tilde {y})} {\rm d}\tilde {y}=-\int _{0}^{\bar {c}} \frac {\tilde {y} \dot {\gamma }}{G_\xi (\bar {c})}{\rm d}\tilde {y}=\frac {(\int _{0}^{\bar {c}} w {\rm d}\tilde {y})=1}{G_\xi (\bar {c})}. \end{equation}
Replacing
$G_\xi (\bar {c})$
in (2.37) with (2.36) leads to
\begin{equation} {I}_1(\bar {c};.)=\frac {1}{G_\xi (\bar {c})}=\sqrt {\frac {\kappa _2^{1/n_2}\left(\dfrac {1}{n_2}+1\right)\left(\dfrac {1}{n_2}+2\right)}{(\bar {c}G_\xi -\tau _{Y,2})_+^{\frac {1}{n_2}+1}(\bar {c}G_\xi \left(\dfrac {1}{n_2}+1\right)+\tau _{Y,2})}}. \end{equation}
To verify this approach, we consider when the yield stress of the displacing fluid is negligible, in which case (2.36) and (2.38) can be simplified further to give
\begin{align} G_\xi =\frac {1}{a\bar {c}^{2n_2+1}+b}, \qquad I_1(\bar {c};.)= (a\bar {c}^{2n_2+1}+d), \end{align}
where
$a=1/{((1/{n_2}+2)\kappa _2^{({1}/{n_2})})^{n_2}}$
and
$d$
is a constant, e.g. for Newtonian fluids,
$G_\xi \approx ({1}/{a\bar {c}^{3}+d})$
. Figure 8(a,b) displays the variation of the axial modified pressure gradient (
$G_\xi$
) and the mean mobility (
$I_1(\bar {c})$
), with concentration for various displacing fluid rheological properties, with the yield stress and the power-law index of the displaced fluid fixed at
$\tau _{Y,1}=0.94$
and
$n_1=0.5$
. The green lines illustrate the curves fitted to expressions of the form
$1/[a\bar {c}^{2n_2+1}+d]$
for the modified pressure gradient, for
$\bar {c}\geqslant \bar {c}_{\textit{min}}$
. We find excellent agreement. Note that
$\bar {c}_{\textit{min}}$
can be calculated analytically with
$\bar {c}_{\textit{min}} = [H/(a\tau _{Y,1})]^{({1}/{2n_2+1})}$
for power-law displacing fluids. This is obtained when
$G_\xi (\bar {c}_{\textit{min}})$
is equal to
$\tau _{Y,1}/H$
, which is required to overcome the yield stress of the displaced fluid. As the rheological properties of the displaced fluid are identical in all the investigated cases,
$G_\xi (\bar {c}_{\textit{min}}) = 0.94$
, at all
$\bar {c}_{\textit{min}}$
values in figure 8(a).
Variation of
$G_\xi (\bar {c})$
(a) and the mean mobility
$I_1(\bar {c})$
(b) by
$\bar {c}$
for isodense flows (
$b=0$
) for invasion of Newtonian, power-law (
$n_2=0.6$
and
$n_2=0.4$
) and Herschel–Bulkley fluids (
$n_2=0.5$
,
$\tau _{Y,2}=0.23$
) into the displaced fluid with a yield stress of
$\tau _{Y,1}=0.94$
and power-law index of
$n_1=0.5$
.
From figures 7(a) and 8(b), we can note that
$I_1\approx a\bar {c}+d$
for
$\bar {c}_{\textit{min}}\leqslant \bar {c}$
when
$\tau _{Y,2}=0.23$
and
$n_2=0.5$
. This is only compatible with the power fluid expressions in the limit of
$n_2 \to 0$
. This gives a hint for where the linear behaviour might come from.
We assume that, for some range of
$\bar {c} \gt \bar {c}_{\textit{min}}$
, particularly when the fluid 1 yield stress is close to 1, meaning that the shear stress at the interface (
$\bar {c}G_\xi$
) is only slightly larger than
$\tau _{Y,2}$
. Using (2.36) we find that
\begin{equation} \bar {c}G_\xi -\tau _{Y,2}\approx \left[\kappa _2^{\frac {1}{n_2}}\left(\frac {1}{n_2}+1\right)\frac {\tau _{Y,2}}{\bar {c}^2}\right]^{\frac {n_2}{n_2+1}} \ll \tau _{Y,2}. \end{equation}
On substituting into (2.38) this can be simplified to
\begin{equation} I_1(\bar {c};.)=\sqrt {\frac {\kappa _2^{\frac {1}{n_2}}\left(\dfrac {1}{n_2}+1\right)\left(\dfrac {1}{n_1}+2\right)}{\kappa _2^\frac {1}{n_2}\left(\dfrac {1}{n_2}+1\right)\dfrac {\tau _{Y,2}}{\bar {c}^2}\left[\tau _{Y,2}\left(\dfrac {1}{n_2}+2\right)+O(\epsilon )\right]}}\approx \frac {\bar {c}}{\tau _{Y,2}}, \end{equation}
i.e. the linear regime observed corresponds to a thin layer of yielded displacing fluid with a static wall layer of displaced fluid.
Figure 9(a,b) illustrates the impact of the power-law index and the yield stress of the displacing fluid on the variation of
$G_\xi$
with concentration. Comparing figure 9(a) with figure 9(b), we observe a better matching between the fitted curve and the calculated
$G_\xi (\bar {c})$
at
$\bar {c}\geqslant \bar {c}_{\textit{min}}$
when the yield stress of the displacing fluid increases (e.g.
$\tau _{Y,2}=0.33,\,0.71$
). The inset of figure 9(a,b) shows the matching between the fitted curves with
$G_\xi$
for
$\bar {c}\geqslant \bar {c}_{\textit{min}}$
. For a small yield stress of the displacing fluid and various power-law indexes, the fitted curves deviate slightly from the calculated
$G_\xi$
close to
$\bar {c}=1$
. However, this deviation is reduced when the yield stress of the displacing fluid is increased (e.g.
$\tau _{Y,2}=0.33,\,0.71$
). The larger yield stresses correspond to cases with a relatively thin layer of yielded fluid 2 close to the interface.
Variation of
$G_\xi$
with
$\bar {c}$
for isodense flows (
$b=0$
) during invasion of Herschel–Bulkley fluids into the displaced fluid with a yield stress of
$\tau _{Y,1}=0.94$
and power-law index of
$n_1=0.5$
. (a) Displacing fluid with (i)
$n_2=0.7$
, (ii)
$n_2=0.5$
, (iii)
$n_2=0.3$
and
$\tau _{Y,2}=0.23$
. (b) Displacing fluid with the same rheological properties as figure 7(a).
3. Planar displacement flows with the D2DGA model
To build our physical intuition on a flow of practical importance, we consider the simplest possible flow: that of a planar displacement. Since the concentrations are always dispersing in the D2DGA model, we must first define what is meant by a planar displacement. The criterion used is that
$\bar {{c}} \approx \bar {{c}}(\xi ,t)$
. For the 2DGA model, this occurs frequently, in so-called steady-state displacements, even in cases of strong inclination and eccentricity. For the D2DGA model, this possibility is less clear. We focus on a scenario where the annulus is both uniform in
$\xi$
and is vertical. We assume that the streamlines are quasi-parallel:
\begin{equation} |\bar {v}| \ll |\bar {w}| \,\,\,\,\Longleftarrow \,\,\,\,\, \left | \frac {\partial \varPsi }{\partial \xi } \right | \ll \left | \frac {1}{r_a}\frac {\partial \varPsi }{\partial \phi } \right |. \end{equation}
For clarity, we shall assume that
$|\bar {v}| = O(\epsilon )$
for some unspecified streamline aspect ratio
$\epsilon \ll 1$
. At leading order we can simplify (2.13) to
\begin{equation} H\bar {w} = H^3[I_1(\bar {c}) G_\xi -I_2(\bar {c}) G_{b,\xi }] = H^3 \left [ - I_1(\bar {c}) \left (\frac {\partial p}{\partial \xi }+\frac {\rho }{{\textit{Fr}}^{*2}} \right ) - b r_a I_2(\bar {c}) \right ]. \end{equation}
The azimuthal buoyancy term is zero and the azimuthal pressure gradient is also
$O(\epsilon )$
to be compatible with
$|\bar {v}| = O(\epsilon )$
.
The invasion of fluid 2 into a Hele-Shaw cell filled with fluid 1 (a). Here N and W represent the narrow and wide sides of the cell and the colour displays the gap-averaged concentration. The evolution of the finger into the dispersive front at the (
$\phi , \xi$
) plane is illustrated in panel (b).
Turning to the transport problem (2.19), we see that (2.21) simplifies to
\begin{equation} (q_\phi ,q_\xi )=r_a H q_0(\bar {c}) (\bar {v},\bar {w}) + r_a b H^3 I_3(\bar {c}) (0,1). \end{equation}
If we write
\begin{equation} \bar {c}(\phi ,\xi ,t)=\bar {c}_0(\xi ,t)+\epsilon \bar {c}_1(\phi ,\xi ,t), \end{equation}
and substitute into (2.19) we find that
\begin{align} 0 & = \frac {\partial [H r_a \bar {c}_0]}{\partial t } + \frac {\partial }{\partial \xi } [r_a H \bar {w} q_0(\bar {c}_0) + r_a b H^3 I_3(\bar {c}_0)] + \frac {\partial }{\partial \phi } [H \bar {v} q_0 (\bar {c}_0)] \nonumber \\ &\quad +\epsilon \frac {\partial }{\partial \xi } [r_a H \bar {w} q_0^\prime (\bar {c}_0)\bar {c}_1 + r_a b H^3 I_3^\prime (\bar {c}_0)\bar {c}_1] + O(\epsilon ^2), \end{align}
with
$^\prime$
denoting the derivative with respect to
$\bar {c}$
. It is interesting to note that since
$\bar {v}= O(\epsilon )$
,
$r_a =1$
and
$H = H(\phi )$
for a uniform annulus, the leading-order equation is
\begin{equation} 0 = \frac {\partial \bar {c}_0}{\partial t } + \frac {\partial }{\partial \xi } [ q_0(\bar {c}_0) + b H^2 I_3(\bar {c}_0)] . \end{equation}
This equation allows for a parametric
$\phi$
dependence, i.e. due to
$H(\phi )$
, while the advective dispersion is only in the
$\xi$
direction. Potentially, this might allow conditions where the dispersive behaviour might be suppressed, e.g. by buoyancy and in a different manner for large and small
$H$
.
A first-order perturbation solution might also be calculable for
$H = H(\phi )$
. The terms involving
$\bar {c}_1$
are also advective and
$\bar {v}$
can be reconstructed from
$\bar {w}$
, using incompressibility. However, instead we proceed assuming that
$H(\phi )$
is near uniform, e.g. if the eccentricity
$e = O(\epsilon )$
:
$H \sim 1$
. The
$\phi$
-derivative term is then
$O(\epsilon ^2)$
. We may consider that the fluids are pumped at speed
$\bar {w} = 1$
(effectively due to the scaling adopted). The leading-order evolution equation for
$\bar {c}_0$
becomes
\begin{equation} 0 = \frac {\partial \bar {c}_0}{\partial t } + \frac {\partial }{\partial \xi } [ q_0(\bar {c}_0) + b I_3(\bar {c}_0)] . \end{equation}
3.1. Characterising the dispersive front
Equation (3.7) is a quasilinear hyperbolic equation. Using kinematic wave theory (Whitham Reference Whitham1974), the front velocity can be calculated by (3.8) provided that
$\bar {c}(\xi ,t)$
remains differentiable:
\begin{equation} w_f(\bar {c}_0)= q_0^\prime (\bar {c}_0) + b I_3^\prime (\bar {c}_0). \end{equation}
We recall that the premise of the D2DGA model is that a central finger of fluid 2 disperses along the channel. Provided
$w_f(\bar {c})$
decreases monotonically with
$\bar {c}$
, the concentration will disperse. In other situations a shock may form, e.g. between two concentrations
$\bar {c}_m$
and
$\bar {c}_M$
. The interval
$[\bar {c}_m,\bar {c}_M]$
moves with constant shock velocity
$w_s$
, which is also the limiting value of
$w_f(\bar {c})$
from each end. The shock velocity and its height (
$\bar {c}_M-\bar {c}_m$
) can be determined by the equal-areas rule or the Rankine–Hugoniot conditions (Whitham Reference Whitham1974). Similar methods have been applied in the cementing context by Pelipenko & Frigaard (Reference Pelipenko and Frigaard2004c
), Zare et al. (Reference Zare, Roustaei and Frigaard2017).
In order to explore the behaviour of solutions to (3.7), we use two methods. First, we calculate the isotropic flux function
$q_0(\bar {c}_0)$
and the buoyancy flux distribution
$I_3(\bar {c}_0)$
by solving for the velocity profile and pressure gradient, imposing the velocity
$(\bar {v},\bar {w}) = (0,1)$
, reflecting the mean imposed speed, as described in Appendix A.2.2. We then analyse these functions to understand the wavespeeds
$w_f(\bar {c})$
and occurrence of shocks. Here, we have benchmarked our solution method against those of Zare et al. (Reference Zare, Roustaei and Frigaard2017), who develop a lubrication theory displacement model for two Bingham fluids displacing along a channel at high Péclet number. Secondly, (3.7) is solved numerically as an evolution problem, here using a total variation diminishing method with a minmod flux limiter to eliminate numerical oscillations at steep concentration gradients (Blunt & Rubin Reference Blunt and Rubin1992; Yang & Yortsos Reference Yang and Yortsos1997).
As examples, figure 11 plots the displacing fluid flux function
$q_0(\bar {c}_0) + b I_3(\bar {c}_0)$
against
$\bar {c}_0$
, and plots
$\bar {c}_0(\xi ,t)$
against the similarity variable
$\xi /t$
. For the latter, once distant from any initial condition,
$\bar {c}_0(\xi ,t)$
collapses quickly to a function of
$\xi /t$
, illustrating both wavespeeds and shocks. Four regimes can be observed for the different rheological properties studied: dispersive front, shock, spike and static wall layers. These regimes are as defined previously for Bingham–Newtonian fluids (Zare et al. Reference Zare, Roustaei and Frigaard2017). In the dispersive regime, the front speed decreases monotonically from the centre of the channel to the walls (no shock). In the shock regime, the front speed is constant at the shock speed and we observe a shock in the interval (
$0\lt \bar {c}_0\lt \bar {c}_M$
), where
$\bar {c}_M$
is the shock height. Similarly, in the spike regime, the front speed is constant between
$\bar {c}_m\lt \bar {c}\lt \bar {c}_M$
, but now
$\bar {c}_m \gt 0$
and the tip of the front (at small
$\bar {c}$
) moves ahead of the shock. In the static wall layer regime, the layer of fluid 1 is not removed (this only occurs for positive yield stress).
Four flow regimes are identified by plotting (a) the displacing fluid flux function
$q_0(\bar {c}_0) + b I_3(\bar {c}_0)$
against
$\bar {c}_0$
, and (b)
$\bar{c}_0$
against the similarity variable (
$\xi/t$
). The dimensionless variables are fixed at
$\bar {w}=1, H=1$
, Dispersion: Newtonian fluids at
$ \boldsymbol{\kappa }=[1\hspace {0.1cm}1],b=0,$
Spike:
$\boldsymbol{\kappa }=[3\hspace {0.1cm}1],\boldsymbol{\unicode{x03C4}_Y}=[0.25\hspace {0.1cm}0],b=50,\boldsymbol{n}=[1\hspace {0.1cm}1],$
Shock:
$ \boldsymbol{\kappa }=[10\hspace {0.1cm}1],\boldsymbol{\unicode{x03C4}_Y}=[12.5\hspace {0.1cm}0],b=100,\boldsymbol{n}=[1\hspace {0.1cm}1],$
Static wall layer:
$\boldsymbol{\kappa }=[1\hspace {0.1cm}1],\boldsymbol{\unicode{x03C4}_Y}=[20\hspace {0.1cm}5],b=0, \boldsymbol{n}=[0.5\hspace {0.1cm}0.5].$
As we move between regimes, we can observe corresponding changes in the flux:
$q_0(\bar {c}_0) + b I_3(\bar {c}_0)$
. The flux function is concave in the dispersive regime. As a shock emerges, it becomes convex over an interval
$[\bar {c}_1,\bar {c}_2]$
, which lies within the shock interval
$[\bar {c}_m, \bar {c}_M ]$
. Equation (3.8) implies that larger concentrations within
$[\bar {c}_1,\bar {c}_2]$
travel faster, eventually leading to the disappearance of the dispersive finger. Thus, to satisfy mass conservation, instead a shock emerges with a velocity equal to the front velocity (3.8) and joining the two points
$\bar {c}_m$
and
$\bar {c}_M$
:
\begin{equation} q(\bar {c}_M)-q(\bar {c}_m)=w_{s}(\bar {c}_M-\bar {c}_m). \end{equation}
As expected, there is good agreement between the shock velocity obtained directly from (3.9) and the functions in figure 11(a), compared with those plotted in figure 11(b), which are computed from the evolution problem.
3.2. Results: dispersive flow regimes
We now explore more generally the parametric dependency of the four identified dispersive flows. We first compute the front speed
$w_f(\bar {c}_0)$
at
$\bar {c}_0=0$
for both Newtonian and shear-thinning fluids for a wide range of buoyancy number (
$b$
), viscosity ratio (
$m$
), and power-law index (
$n_k$
). Our main goal here is to find a critical
$b$
value to reduce the impact of dispersion and accordingly result in a more stable front in the range studied.
Figure 12 explores the impact of
$m=\kappa _1/\kappa _2$
and
$b$
on the front speed for two Newtonian fluids. The shaded contours give the centreline front speed
$w_f(\bar {c}_0=0)$
. Note that the effect of the buoyancy number and the viscosity ratio on
$q_0^{\prime}(\bar {c}_0)+ b I_3^{\prime}(\bar {c}_0)$
is negligible at
$\bar {c}_0=0$
, and it maintains a constant value of 1.5 at the centre of the channel for Newtonian fluids, coming from the maximum velocity of the Poiseuille flow between two parallel plates. However, when shock develops, the shock velocity might be greater (or less) than
$q_0^{\prime}(0)+ b I_3^{\prime}(0)$
and
$w_f$
can be changed by both
$b$
and
$m$
.
Shaded contours of the front speed (
$w_f(\bar {c}_0)$
) at
$\bar {c}_0=0$
for Newtonian fluids with (a)
$b=0$
, (b)
$b=10$
, (c)
$b=100$
, (d)
$b=1000$
. The contour lines in (a) represent the front speed and the broken red lines correspond to the critical viscosity ratios of Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999). The subfigures (middle and right-hand side plots) give examples of
$q_0^\prime (\bar {c}_0) + b I_3^\prime (\bar {c}_0)$
and consequent
$w_f(\bar {c}_0)$
, at indicated values of
$(\kappa _1,\kappa _2)$
. The first column of subfigures (middle plots) explore
$m \geqslant 1$
(viscosity-unstable flow) and the second column of subfigures (right-hand side plots) explore
$m \leqslant 1$
(viscosity-stable flow).
Generally speaking, the shock velocity approaches the mean velocity (
$w_f(\bar {c}_0) \to 1$
), and the front becomes more stable when
$b \gg 1$
; see figure 12. There is less dispersion for larger
$b$
and larger
$m$
. Figure 12(a) shows the transition from dispersive front to shock regime based on the viscosity ratio of fluids at
$b=0$
. This transition happens at a critical viscosity ratio (
$M^{\textit{min}}_3=1.5$
), which is identical with the predictions of Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999), who considered Newtonian fluids only. To translate from our dimensionless parameters to the
$(U,M)$
of Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999), we have
$b=3m^{0.5}/U$
and
$m=M$
. Here
$M^{\textit{max}}_1$
and
$M^{\textit{min}}_3$
represent the upper limit for the dispersion regime and the lower limit for the shock regime, respectively. In the region between
$M^{\textit{min}}_3$
and
$M^{\textit{max}}_1=2.25$
, the shock speed is equal to
$1.5$
and exceeds 1.5 when
$m=M\gt M^{\textit{max}}_1$
.
Examples of the effect of the viscosity ratio on the front speed are illustrated as subfigures in figure 12 for Newtonian fluids. In the second subfigure of figure 12(a), a more viscous fluid is injected into a less viscous fluid; hence,
$m\leqslant 1$
and flow is viscosity stable. As the viscosity ratio reduces, the front velocity at large concentration increases, implying that the wall’s residual fluid starts to move at a higher speed. In the first subfigure of figure 12(a), the viscosity ratio between the displacing and the displaced fluid is unfavourable, meaning that the displacing fluid is less viscous than the displaced fluid. The front is dispersive at
$m=1$
and a shock develops at
$m=2,\,3,\,4,\,5$
.
The increase in the front speed results in a more unstable regime, making the primary cementing operations less effective. In figure 12(b–d) the displacing fluid is denser than the displaced fluid, i.e.
$b\gt 0$
. Similar to the defined regions of the (
$U,M)$
plane in (Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999), at low
$U$
(or high
$b$
) values, the flow regime is mainly spike. In figures 12(c) and 12(d), buoyancy completely suppresses the spike regime as
$b$
increases from
$100$
to
$1000$
.
Figure 13 shows the front speed at
$\bar {c}_0=0$
for two shear-thinning fluids at different
$n_k$
and
$b$
. We have fixed
$\kappa _1 = \kappa _2 = 1$
. Similar to the Newtonian fluids, in figure 13 we investigate the influence of the buoyancy number on the front speed for shear-thinning fluids. The effects are somewhat similar to those for Newtonian fluids. The general trend is that increasing the buoyancy number can improve the displacement process as the front velocity at
$\bar {c}_0=0$
approaches the mean velocity. This also happens as
$n_1$
is reduced, simply because the displaced fluid becomes progressively shear thinning with a plug-like velocity profile.
Front speed (
$w_f(\bar {c}_0)$
) at
$\bar {c}_0=0$
for power-law fluids with (a)
$b=0$
, (b)
$b=10$
, (c)
$b=20$
, (d)
$b=70$
. The shaded contour lines represent the front speed. The subfigures (middle and right-hand side plots) give examples of
$q_0^\prime (\bar {c}_0) + b I_3^\prime (\bar {c}_0)$
and consequent
$w_f(\bar {c}_0)$
, at indicated values of
$(n_1,n_2)$
. The first (middle plots) and second (right-hand side plots) column of subfigures explore the impact of the power-law index of displacing and displaced fluids on the front speed, respectively.
For
$b=0$
, the main transition in flow type is between dispersive and spike, as
$n_2$
decreases; see figure 13(a). At larger
$b$
the dispersive regimes are suppressed, there are some spike regime flows, but these too are suppressed in favour of shock regimes at larger
$b$
; see figure 13(b–d).
These trends at
$b=0$
are illustrated more clearly in the subfigures of figure 13(a). For a fixed
$n_2$
, the front height increases (or the front speed reduces) as
$n_1$
decreases; see the second subfigure of figure 13(a). This trend has been previously observed, where the front velocity reduction is attributed to the downstream plug velocity, which can limit the velocity of the displacement front (Taghavi et al. Reference Taghavi, Seon, Martinez and Frigaard2009; Zare et al. Reference Zare, Roustaei and Frigaard2017). Conversely, the front velocity can increase, and flow can become less stable as
$n_2$
decreases for a fixed
$n_1$
; see the first subfigure of figure 13(a). This clearly stems from the lower effective viscosity of the displacing fluid compared with the displaced fluid.
3.3. Stability of the dispersive front: Muskat analysis
We now explore the stability of the dispersive planar front to an advancing narrow finger of fluid 2 in the
$(\phi ,\xi )$
plane. Our dispersive model is built on the premise that a dispersive finger advances along the centre of the annular gap, so that a finger in the
$(\phi ,\xi )$
plane appears slightly contradictory. Our working assumption is that (3.7) describes the base flow completely, which leads to the dispersing planar flow. We then simply question the stability of this mathematical system. The underlying approach is attributed to Muskat (Reference Muskat1938). A similar approach has been applied by Pascal (Reference Pascal1984) in considering stability of a frontal porous media displacement and also by Pelipenko & Frigaard (Reference Pelipenko and Frigaard2004c
) in analysing stability of the 2DGA model. However, neither has considered a dispersive front.
The uniform velocity assumption is applied to (3.2), which means that the underlying pressure gradient as we move through the dispersing front is controlled by
$\bar {c}_0$
, and given by
\begin{equation} \frac {\partial p}{\partial \xi } = -\frac {\rho (\bar {c}_0)}{{\textit{Fr}}^{*2}} - \frac {1 + b I_2(\bar {c}_0)}{I_1(\bar {c}_0)}. \end{equation}
Assuming that the finger is thin in the
$\phi$
direction, the pressure transmits across the finger:
\begin{equation} \frac {\partial p_{\textit{finger}}}{\partial \xi } = -\frac {\rho (\bar {c}_0)}{{\textit{Fr}}^{*2}} - \frac {1 + b I_2(\bar {c}_0)}{I_1(\bar {c}_0)}. \end{equation}
However, the finger is full of fluid 2. Hence, the velocity of the fluid within the finger
$\bar {w}_{\textit{finger}}$
is given by (3.2) with
$\bar {c}_0=1$
and the pressure gradient from (3.11):
\begin{eqnarray} \bar {w}_{\textit{finger}}(\bar {c}_0) & = & = - I_1(1) \left (\frac {\partial p_{\textit{finger}}}{\partial \xi }+\frac {\rho (1)}{{\textit{Fr}}^{*2}} \right ) - b I_2(1) \nonumber \\ &=& \frac {I_1(1)}{I_1(\bar {c}_0)} + b I_1(1) \left ( \frac {I_2(\bar {c}_0)}{I_1(\bar {c}_0)} + {\bar {c}_0}-1\right ). \end{eqnarray}
Note that we have used the fact that
$I_2(1)=0$
in the above. We observe that the finger velocity is driven both by a form of mobility ratio
$I_1(1)/I_1(\bar {c}_0)$
(first term) and by buoyancy (second term).
Having found the finger velocity, the usual approach is to compare with a suitable wavespeed to determine stability or otherwise. However, in our case the front is dispersing with speed
$\bar {w}_{f}(\bar {c}_0)$
. We consider that
$\bar {w}_{f}(\bar {c}_0)$
has been computed to include any segments of
$\bar {c}_0$
over which there is a shock front, as described in § 3.1. We then define
\begin{equation} \Delta w(\bar {c}_0) = \bar {w}_{\textit{finger}}(\bar {c}_0) - \bar {w}_{f}(\bar {c}_0) \gt 0. \end{equation}
The simplest regime to classify is when the flow is unstable. For an unstable regime,
$\triangle w(\overline{c}_{0} ) \gt0 $
for all
$\overline{c}_0\in [0, 1] $
. In this case, the finger of fluid 2 advances through the dispersive front at all concentrations.
The notion of stability is harder to define. Conceptually at least, we can conceive that the finger might not penetrate the dispersing front at all or might penetrate for some range of concentrations. Loosely, we call these eventualities stable and partial penetration. These are illustrated schematically in figure 14(a).
(a) Schematic of flow classification. (b) Schematic comparison of the finger velocity with the front velocity for Newtonian fluids on the (
$b,m$
) plane at three distinct flow regimes: (i) stable (green), (ii) partial penetration (blue) and (iii) unstable (red) regimes. (c) The frontal speed (
$w_f(\bar {c})$
) is decoupled into two parts: (i)
$q_0^{\prime}(\bar {c}_0)$
and
$b*I_3^{\prime}(\bar {c}_0)$
and then plotted on the (
$b,m$
) plane.
To aid in classification, first note that
$w(\bar {c}_0)$
, after processing to take account of shocks, is always a decreasing function of
$\bar {c}_0$
. On the other hand, little can be said in general regarding the behaviour of
$\bar {w}_{\textit{finger}}(\bar {c}_0)$
. There may be one or more concentrations for which
$\Delta w(\bar {c}_0) = 0$
. We define
$c_{\textit{critical}}$
to be the largest value of
$\bar {c}_0 \in [0,1]$
for which
$\Delta w(\bar {c}_0) = 0$
. Observe that, as
$\bar {c}_0 \to 1$
, the derivatives of both flux functions
$q_0$
and
$I_3$
vanish, so that
$w_f(1) = 0$
. On the other hand, since as
$\bar {c}_0 \to 1$
,
$I_2(\bar {c}_0) \to 0$
, we can see from (3.12) that
$\bar {w}_{\textit{finger}}(1) = 1$
. Therefore,
$\Delta w(1) = 1$
and the finger will always start to penetrate the dispersive front. If there is no
$c_{\textit{critical}}$
then it follows that
$\Delta w(\bar {c}_0) \gt 0$
for all
$\bar {c}_0 \in [0,1]$
and the flow is unstable.
In terms of specifying a stable flow, a necessary condition is that there exists a
$c_{\textit{critical}} \in [0,1)$
. One might select a threshold concentration such that, e.g.
$c_{\textit{critical}} \geqslant 0.9$
designates stability, i.e. the finger is allowed to penetrate through at most 10 % of the dispersive front. In this way the influence of buoyancy in stabilising the flow can be assessed. Increasing
$b$
tends to suppress
$\bar {w}_{\textit{finger}}$
. Also, for
$\bar {c}_0$
close to
$1$
, it is found that
$I_3^\prime (\bar {c}_0) \gt 0$
. Thus, although
$q_0^\prime (1) = 0$
, sufficiently large
$b$
results in
$w_f(\bar {c}_0) \gt 0$
in an interval close to
$1$
, i.e. large enough
$b\gt 0$
can result in
$\Delta w(\bar {c}_0) \lt 0$
and stability for any fixed
$\bar {c}_0 \lt 1$
. This motivates the idea of adopting a threshold value, but the actual threshold adopted lacks meaning.
Therefore, we turn to the mean displacement speed in order to identify a threshold. We denote by
$\bar {c}_*$
the value of
$\bar {c}_0$
for which
$w_f(\bar {c}_0) = 1$
, i.e.
$\bar {c}_*$
is the concentration that disperses at the mean speed. Note that
$w_f(\bar {c}_0) - 1$
integrates to zero over
$[0,1]$
, i.e.
$\bar {c}_*$
is a volume-weighted mean concentration: the same volume disperses ahead of
$\bar {c}_*$
as behind
$\bar {c}_*$
. Indeed if there is shock moving at speed
$1$
, we would need to define
$\bar {c}_*$
in this way. With the above explanation and meaning, we define stable and partial penetration as follows.
For a stable regime, the displacement is regarded as stable if
$\Delta w(\bar {c}_0)\lt 0$
for all
$\bar {c}_0 \in [0,c_{\textit{critical}}]$
and
$c_*\lt c_{\textit{critical}}$
. Thus, fingers of fluid 2 do not penetrate ahead of the mean pumping speed, which is dispersing at concentration
$\bar {c}_*$
. Fingers will penetrate in the interval
$(c_{\textit{critical}},1]$
, but at least 50 % of the front is unaffected. As the front advances, the penetrated interval
$(c_{\textit{critical}},1]$
will recede behind the front advancing at the mean speed.
For a partial penetration regime,
$\Delta w(\bar {c}_0)\gt 0$
for
$\bar {c}_0 \gt c_{\textit{critical}}$
and
$c_*\geqslant c_{\textit{critical}}$
. In other words, the finger penetrates through at least 50 % of the front, but fails to penetrate through to the pure fluid 1 layer, halting at some critical value.
Flow regimes classification based on Muskat’s analysis for Newtonian fluid flows at various viscosity ratios (
$m$
) and buoyancy numbers (
$b$
): (a)
$b=0$
, (b)
$b=10$
, (c)
$b=100$
, (d)
$b=1000$
. Red triangles denote unstable; blue squares denote partial penetration; green circles are stable.
Flow regimes classification based on Muskat’s analysis for power-law fluid flows at various power-law index (
$n_k$
) and buoyancy numbers (
$b$
): (a)
$b=0$
, (b)
$b=10$
, (c)
$b=20$
, (d)
$b=70$
. Red triangles denote unstable; blue squares denote partial penetration; green circles are stable.
(a,c) Front speed (
$w_f(\bar {c}_0)$
) is plotted against
$\bar {c}_0$
at
$b=1,\,\kappa _1=0.3,\, \kappa _2=0.1$
. The colour map represents
$\Delta w(\bar {c}_0)$
for (b) various
$\bar {c}_0$
and
$\tau _{Y,1}$
at fixed
$\tau _{Y,2}=0$
and for (d) different
$\bar {c}_0$
and
$\tau _{Y,2}$
at fixed
$\tau _{Y,1}=0.05$
. The inset figures give examples of a comparison between finger and front speeds.
3.4. Results
For Newtonian fluid pairs, due to linear superposition buoyancy effects are not closely coupled to viscous effects. Figure 14(b) presents variations in stability classification due at nine different
$m$
and
$b$
values. Broadly speaking, small
$m$
and larger
$b$
enhance stability. The limit values and monotonicity discussed above are illustrated in these subfigures, as is the treatment of the
$q_0^{\prime}(\bar {c}_0)+ b I_3^{\prime}(\bar {c}_0)$
in order to determine shocks. In figure 14(c) we plot the individual contributions
$q_0^{\prime}(\bar {c}_0)$
and
$b I_3^{\prime}(\bar {c}_0)$
. From this we can see how buoyancy comes to dominate
$q_0^{\prime}$
, increasing
$w_f$
.
Now, turning to the results of the Muskat analysis/classification of the flow regimes, we compared the front velocity with the finger speed for all ranges of concentrations. Then, we classified them based on the explained criteria. Figure 15 displays the flow classification based on Muskat’s analysis for two Newtonian fluids. In the isodense displacements, the front becomes unstable for values of
$m$
higher than a cutoff
$m_c$
, where
$m_c\sim 1.5$
. This implies that the finger velocity is higher than the front speed at
$\bar {c}=0$
when
$m\geqslant m_c$
, as shown in figure 15(a). Under a favourable density ratio between two fluids, i.e.
$b\gt 0$
, the unstable region reduces, and finally, at
$b=1000$
, the buoyancy force dominates, and the front becomes stable for all ranges of the viscosity ratio.
At high values of the viscosity ratios, a much larger buoyancy number is required to stabilise the flow. Therefore, a transition from unstable regime to stable occurs, and the mud removal process becomes more effective when the displacing fluid becomes denser and more viscous than the displaced fluid, as is indeed physically intuitive.
The influence of the power-law index values on the frontal behaviour of two shear-thinning fluids is illustrated in figure 16. Essentially, as the shear thinning of the displaced fluid decreases and the shear-thinning fluid of the displacing fluid increases, i.e.
$n_1\to 1$
and
$n_2\to 0.2$
, the front becomes less stable and the displacement effectiveness reduces. Furthermore, as the effective viscosity of fluids is smaller than in the previous scenario (Newtonian fluids), the critical buoyancy number required to stabilise the flow (
$b \approx 70$
) is far below the threshold for Newtonian flows (
$b\approx 1000$
).
Finally, we explore some of the effects of yield stress on the displacement. Including a yield stress creates more possibilities for either fluid to be unyielded over part or all of the layer thickness, which can combine with buoyancy in either direction to produce many effects. Therefore, the results presented are primarily intended as examples.
In figure 17 we fix a modest stabilising buoyancy
$b=1$
, set
$n_1=n_2 = 1$
and consider an unstable ratio of the consistencies. Figures 17(a) and 17(b) show the results of increasing
$\tau _{Y,1}$
from zero, keeping
$\tau _{Y,2} = 0$
. Initially the displacement has a spike flow, but as
$\tau _{Y,1}$
increases the spike changes to a shock. The velocity profile develops a static wall layer of fluid 1. The shock velocity increases with
$\tau _{Y,1}$
in this case (panel a). Considering now the effects on stability, we have plotted
$\Delta w(\bar {c}_0)$
against
$\tau _{Y,1}$
in the form of a colour map. We see that initially for low
$\tau _{Y,1}$
we have a band of low concentrations for which
$\Delta w(\bar {c}_0) \lt 0$
, but this does not extend to
$\bar {c}_*$
. Consequently, these flows show partial penetration and are not stable. As the yield stress is increased, we find that
$\Delta w(\bar {c}_0) \gt 0$
for all
$\bar {c}_0$
and the flows become unstable. This is illustrated within the inset figures of figure 17(b),
In figures 17(c) and 17(d) we explore instead the effect of increasing
$\tau _{Y,2}$
from zero, keeping
$\tau _{Y,1} = 0.05$
. In terms of flow type (panel c), we always find spike-type frontal displacements, with the spike region growing larger as
$\tau _{Y,2}$
increases. As the shock front speed moves below
$1$
there is a sudden change in
$\bar {c}_*$
, which drops below
$0.5$
. At the same time,
$\bar {c}_{\textit{critical}}$
steadily increases so that we transition from partial penetration to a stable fingering regime (panel d).
4. Conclusions
The main focus of this study has been the extension of the D2DGA approach of Zhang & Frigaard (Reference Zhang and Frigaard2022) to a wide range of non-Newtonian fluids, such as are commonly used in primary cementing flows. The underlying approach is illustrated in figure 3, namely developing a closure expression for the functional relationship of the pressure gradient to the gap-averaged velocity (streamfunction), based on the assumption that the fluids are layered across the annular gap. As shown in Zhang & Frigaard (Reference Zhang and Frigaard2022, Reference Zhang and Frigaard2023), this approach leads to a greatly improved physical representation of dispersion that takes place on the scale of the annular gap. Although there is no reason to suspect that validity of the D2DGA approach fails due to adoption of different generalised Newtonian fluid rheologies, we admit that external verification is an important missing step so far. Computation of three-dimensional annular displacement flows using Herschel–Bulkley fluid models is underway, with the purpose of making these comparisons. We also have ongoing experimental displacement work that may eventually be used for this purpose. This will complete the modelling framework developed here.
We have developed the underlying closure expressions needed to compute the streamfunction, given in terms of the mobility functions
$I_1$
and
$I_2$
, and explored their parametric variations. Unlike the Newtonian fluid pairs, an analytical solution to find the velocity field is no longer possible. Instead it is necessary to compute the velocity numerically, and here we have used an augmented Lagrangian method. The velocity solutions can have plug regions of unyielded fluid, either in the centre of the annular gap or at the walls. Also, the shear stress direction changes across the gap width, due to the different effects of buoyancy in each layer. In general, the result of finding the closure expressions is that the mean mobility
$I_1 \gt 0$
and (2.16) remains elliptic. The solution method for the gap-scale velocity is detailed in Appendix A.
We have also derived the general expressions needed for evolution of the displacing fluid concentration
$\bar {c}$
. The volumetric flux function splits into two parts: a flux that partitions the total areal flow rate locally and is in the direction of the mean velocity (
$\boldsymbol{q}_0$
), and a buoyancy driven term
$bI_3$
. This differs from the Newtonian case in that the two components are coupled through the streamfunction as well as the concentration, i.e. the linearity is lost. The resulting nonlinear evolution problem closure expressions must be evaluated computationally from the local gap-scale velocity solution.
The need to calculate each closure expression numerically makes computation of the complete D2DGA displacement flow impractical without, e.g. tabulating each function. Instead, to gain physical insight into the effects of dispersion, we have analysed a planar displacement flow using the D2DGA model. The analysis consists of two parts. First, we are able to compute the dispersive wavespeeds of the planar displacement front. In the Newtonian case, this incorporates the earlier results of Yang & Yortsos (Reference Yang and Yortsos1997), Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999), but builds on this for power-law and Herschel–Bulkley fluids. Similar to the displacement flow of Bingham fluids by Newtonian fluids (Zare et al. Reference Zare, Roustaei and Frigaard2017), frontal behaviours in Herschel–Bulkley fluid displacements can be categorised into (a) dispersion, (b) spike, (c) shock and (d) static wall layers. Fluid rheological properties and buoyancy control the transition from one regime to another. Shear thinning in the displaced fluid tends to reduce dispersion, while shear thinning of the displacing fluid has little effect. Increasing the yield stress of the displacing fluid enhances the velocity of the displaced fluid layers, without necessarily changing the flow regime. As the yield stress of the displaced fluid is increased, eventually the static wall layer regime arises. The effects on dispersion are less clear.
Including positive buoyancy
$b\gt 0$
tends to reduce dispersion, moving from fully dispersive flows, through spike regimes to central shocks of widening extent. The suppression of dispersion and the shock development, on applying a positive buoyancy, stems from the change of the flux function (
$q_0(\bar {c})+ b I_3(\bar {c})$
) from concave into convex at small values of
$\bar {c}$
. The height and velocity of the shock are determined using the Rankine–Hugoniot conditions. On the other hand, shock speeds are always positive regardless of how large
$b$
is, due to the behaviour of the buoyancy flux distribution
$I_3$
, as
$\bar {c}\to 0$
and
$1$
.
Secondly, we conducted a dispersive Muskat-style analysis to calculate the velocity of a narrow finger penetrating the dispersing front. Since in the Muskat analysis it is assumed that the finger is very thin, the pressure gradient within the finger is the same as that in the external fluid, and accordingly, the finger velocity depends on
$\bar {c}$
. We classify the behaviour as unstable, when the finger penetrates through the front and into the layer of fluid 1. The finger velocity is always faster than the front velocity at
$\bar {c}=1$
, so the finger can always penetrate into the frontal regime. We then find the concentration
$\bar {c}_*$
that disperses at the mean speed. If the finger velocity at
$\bar {c}_*$
does not exceed the mean velocity, we regard the flow as stable. Otherwise, the flow is classified as a partial penetration.
Applying this analysis to frontal displacement flows of two Newtonian fluids gives results that are somewhat intuitive. Increasing displacing fluid viscosity tends to stabilise the front. For an isodense displacement, unstable flows are found for
$m\gt 1.5$
. Adding buoyancy to the displacement (
$b\gt 0$
) progressively stabilises the displacement. This stabilisation via positive buoyancy is attributed to both suppression of the frontal dispersion and reduction of the finger velocity, by at most
$bI_1(1)$
.
For two power-law fluids with identical
$\kappa _j$
, on varying
$n_1$
and
$n_2$
, we find that reduced
$n_2$
at fixed
$n_1$
can lead to instability. In this case the displacing fluid in the finger feels stress gradients that are controlled by the less responsive fluid 1: it may shear thin and move faster in the centre of the annular gap. Again, as we increase
$b\gt 0$
, the flows become progressively stabilised. The shear-thinning scenario appears to be more sensitive to
$b$
than two Newtonian fluids.
Our future work is directed towards developing a two-dimensional numerical implementation of the non-Newtonian D2DGA model derived here, in such a way that it can be used for systematic computations of cementing displacement flows on the scale of the well. In order to compute the full two-dimensional model, we need to solve (2.16) at each time step and evolve the fluid concentration via (2.19). The majority of this paper has been concerned with formulating the underlying D2DGA model, to show that it has a similar two-dimensional form to the earlier 2DGA models, exploring the closure relationships via the unctions
$I_1$
,
$I_2$
,
$I_3$
and
$q_0$
, and studying the planar displacement flow. We have also studied the underlying
$y$
-dependent two-dimensional shear flow; (see Appendix A), needed in order to evaluate
$I_1$
,
$I_2$
,
$I_3$
and
$q_0$
. Since the computation of these functions is itself time consuming, when eventually solving (2.16) and (2.19), it is likely that these functions need to be evaluated and tabulated so that a more speedy interpolation of their values can be used. Even given this, we note that, for the 2DGA model, (2.16) is properly posed as a convex minimisation problem or the equivalent elliptic variational inequality. It is likely that a similar approach will need developing here for the fully two-dimensional annular flow.
Funding
The authors gratefully acknowledge the financial support of NSERC, Sanjel Energy Services and Geonomic Ltd. via Alliance Project No. ALLRP 572089-22 ‘Well Integrity Canada 2030’.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Solving the two fluid problem
We now describe the solution of the two fluid layered flow in detail. Whereas the model has been derived using variables that are scaled globally, i.e. using process variables relevant to the annular displacement flow, here we consider only the local flow in a plane channel of width
$2H$
. This is governed by the momentum (2.5), (2.6) and the constitutive laws (2.31). In order to reduce the dimension of the parameter space, we define
$\tilde {y} = y/H$
, so that the interfaces are now at
$\tilde {y} = \pm \bar {c}$
. The governing equations are then identical with (2.5), (2.6) and (2.31), except that the parameters and variables are accented with a tilde and
$H = 1$
, i.e. we solve
\begin{eqnarray} 0 &=& \frac {{\rm d}\boldsymbol{\tilde {\unicode{x03C4}}_2}}{{\rm d} \tilde {y}} + \boldsymbol{\tilde {G}} - (1-\bar {c})\boldsymbol{\tilde {G}_b},\,\,\,\,\,\,\tilde {y} \in [0, \bar {c}), \end{eqnarray}
\begin{eqnarray} 0 &=& \frac {{\rm d}\boldsymbol{\tilde {\unicode{x03C4}}_1}}{{\rm d}\tilde {y}} + \boldsymbol{\tilde {G}} +\bar {c}\boldsymbol{\tilde {G}_b},\,\,\,\,\,\,\tilde {y} \in ( \bar {c},1]. \end{eqnarray}
The constitutive laws for fluid
$k$
are
\begin{equation} \begin{array}{lcl} \boldsymbol{\tilde {\unicode{x03C4}}_k} = \displaystyle {\left ( \tilde {\kappa }_k \left |\frac {{\rm d}\boldsymbol{\tilde {u}}}{{\rm d}\tilde {y}}\right |^{n_k-1} + \tilde {\unicode{x03C4}}_{Y,k} \left |\frac {{\rm d}\boldsymbol{\tilde {u}}}{{\rm d}\tilde {y}}\right |^{-1} \right ) \frac {{\rm d}\boldsymbol{\tilde {u}}}{{\rm d}\tilde {y}}}, &\,\,\,\,& \Leftrightarrow |\boldsymbol{\tilde {\unicode{x03C4}}_k}| \gt \tilde {\unicode{x03C4}}_{Y,k}, \\[1ex] \displaystyle {\frac {{\rm d}\boldsymbol{\tilde {u}}}{{\rm d}\tilde {y}}} = 0,& & \Leftrightarrow |\boldsymbol{\tilde {\unicode{x03C4}}_k}| \leqslant \tilde {\unicode{x03C4}}_{Y,k}. \end{array} \end{equation}
At
$\tilde {y}= 0$
, we have
$\boldsymbol{\tilde {\unicode{x03C4}}_2} = 0$
. At
$\tilde {y}= 1$
, we have
$\boldsymbol{\tilde {u}} = 0$
. The stresses and velocity are continuous at
$\tilde {y}= \bar {c}$
. The variables with the tilde accent are identical to the unaccented variables earlier, with the following exceptions:
\begin{equation} \tilde {\kappa }_k = \kappa _k/H^{n_k}, \qquad\boldsymbol{\tilde {G}} = H \boldsymbol{G}, \qquad\boldsymbol{\tilde {G}_b} = H \boldsymbol{G_b}. \end{equation}
Thus, the velocity, yield stresses and power-law indices are unchanged.
A.1. Variational formulation
Aside from special cases (e.g. Newtonian), it is not possible to find an analytical solution. Due to the yield stress, it is common to use a variational formulation. We take a test space
$V$
:
\begin{align} V = \{\boldsymbol{\tilde {v}}\colon \boldsymbol{\tilde {v}}(1) = 0, \boldsymbol{\tilde {v}} \mbox{ continuous at } \tilde {y}=\bar {c} \}. \end{align}
We may expect that
$V$
is a subspace of the Sobolev space
$[W^{1,p}(0,1)]^2$
, where
$p = 1 +\min \{n_1,n_2\}$
, but we proceed only formally here. On taking the dot product of the momentum equations with
$\boldsymbol{\tilde {v}}$
, integrating across
$[0,1]$
and summing, we arrive at
\begin{align} 0 = \left\langle \frac {{\rm d}\boldsymbol{\tilde {\unicode{x03C4}}} }{{\rm d}y} \boldsymbol{\cdot }\boldsymbol{\tilde {v}} \right\rangle + \langle \boldsymbol{\tilde {G}}.\boldsymbol{\tilde {v}} \rangle -(1-\bar {c}) \langle \boldsymbol{\tilde {G}_b}\boldsymbol{\cdot }\boldsymbol{\tilde {v}} \rangle _2 + \bar {c}\langle \boldsymbol{\tilde {G}_b}\boldsymbol{\cdot }\boldsymbol{\tilde {v}} \rangle _1, \end{align}
where the angle bracket notation is as follows:
\begin{equation} \langle f \rangle _2 =\int _{0}^{\bar {c}} f \, {\rm d}\tilde {y}, \qquad \langle f \rangle _1 =\int _{\bar {c}}^1 f \, {\rm d}\tilde {y}, \qquad \langle f \rangle =\int _{0}^{1} f \, {\rm d}\tilde {y} = \langle f \rangle _2 + \langle f \rangle _1 . \end{equation}
The shorthand
$\boldsymbol{\tilde {\unicode{x03C4}}}$
is taken to mean
$\boldsymbol{\tilde {\unicode{x03C4}}_k}$
when in fluid
$k$
.
On integrating by parts, using the boundary and interfacial conditions, and noting that both
$\boldsymbol{\tilde {G}}$
and
$\boldsymbol{\tilde {G}_b}$
are constant, we arrive at
\begin{align} \left\langle \frac {{\rm d}{\boldsymbol{\tilde {v}}}}{{\rm d}\tilde {y}} \boldsymbol{\cdot }\boldsymbol{\tilde {\unicode{x03C4}}} \right\rangle = \boldsymbol{\tilde {G}} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle -(1-\bar {c}) \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle _2 +\bar {c} \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle _1 . \end{align}
To treat the left-hand side, we define the following functionals for fluid
$k$
:
\begin{align}& a_k(\boldsymbol{\tilde {u}}, \boldsymbol{\tilde {v}})= \tilde {\kappa }_k \left\langle \left | \frac {{\rm d} \boldsymbol{\tilde {u}}}{{\rm d}\tilde {y}} \right |^{n_k-1} \frac {{\rm d} \boldsymbol{\tilde {u}}}{{\rm d}\tilde {y}} \boldsymbol{\cdot }\frac {{\rm d} \boldsymbol{\tilde {v}}}{{\rm d}\tilde {y}} \right\rangle _k,\,\,\,\,k=1,2, \end{align}
\begin{align}&\qquad j_k(\boldsymbol{\tilde {v}})= \tilde {\tau }_{Y,k} \left\langle \left | \frac {{\rm d} \boldsymbol{\tilde {v}}}{{\rm d}\tilde {y}} \right |^{n_k-1} \right\rangle _k,\,\,\,\,k=1,2. \\[9pt] \nonumber \end{align}
Then
$a(\boldsymbol{\tilde {u}}, \boldsymbol{\tilde {v}})=a_1(\boldsymbol{\tilde {u}}, \boldsymbol{\tilde {v}})+a_2(\boldsymbol{\tilde {u}}, \boldsymbol{\tilde {v}})$
;
$j(\boldsymbol{\tilde {v}})=j_1(\boldsymbol{\tilde {v}})+j_2(\boldsymbol{\tilde {v}})$
. Using (A9) and (A10), (A8) leads to the following variational inequality:
\begin{eqnarray} a(\boldsymbol{\tilde {u}}, \boldsymbol{\tilde {v}}-\boldsymbol{\tilde {u}}) + j(\boldsymbol{\tilde {v}})-j(\boldsymbol{\tilde {u}}) & \geqslant & \boldsymbol{\tilde {G}} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}}-\boldsymbol{\tilde {u}} \rangle -(1-\bar {c}) \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} - \boldsymbol{\tilde {u}} \rangle _2 +\bar {c} \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} - \boldsymbol{\tilde {u}} \rangle _1 , \nonumber \\ & & \,\,\,\,\,\,\,\, \boldsymbol{\tilde {u}} \in V\quad \forall \boldsymbol{\tilde {v}} \in V. \end{eqnarray}
This second-order elliptic variational inequality is of classical form and can be used to ensure the existence and uniqueness of the solution. Here we simply assume this. The variational inequality (A11) is also equivalent to the following minimisation:
\begin{equation} \underset {\boldsymbol{\tilde {v}}\in V}{\min } \,J(\boldsymbol{\tilde {v}}) = \frac {a_1(\boldsymbol{\tilde {v}},\boldsymbol{\tilde {v}})}{n_1+1} + \frac {a_2(\boldsymbol{\tilde {v}},\boldsymbol{\tilde {v}})}{n_2+1} + j(\boldsymbol{\tilde {v}}) - \boldsymbol{\tilde {G}} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle + (1-\bar {c}) \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle _2 - \bar {c} \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle _1 . \end{equation}
In both variational inequality and minimisation formulations, we observe that the buoyancy vector acts in the opposite directions within each fluid layer.
In practice, we solve the minimisation formulation, using an augmented Lagrangian approach, which we outline in Appendix A.2. Our main aim in finding the solution is so that we define the closure expressions needed for the D2DGA model earlier. The pointwise values of
$\boldsymbol{\tilde {u}}$
are needed to evaluate functions such as
$I_1$
,
$I_2$
,
$q_0$
, etc. There is also direct interest in the mean velocity
$\langle \boldsymbol{\tilde {u}} \rangle$
, which gives
${\nabla }_a \varPsi = 2H (\langle \tilde {u}_{\xi } \rangle , -\langle \tilde {u}_{\phi } \rangle$
.
There are generally two scenarios in which we might wish to compute a solution as part of the D2DGA model. Firstly, it can be that the modified pressure gradient (
$\boldsymbol{\tilde {G}}$
) is specified and we wish to find the velocity
$\boldsymbol{\tilde {u}}$
(hence, the mean velocity and other functionals). The second way in which we need to compute the solution is when
$\langle \boldsymbol{\tilde {u}} \rangle$
is specified, i.e. via
${\nabla }_a \varPsi$
from the D2DGA model. In this case, the modified pressure gradient
$\boldsymbol{\tilde {G}}$
is unknown and must be found along with the solution
$\boldsymbol{\tilde {u}}$
.
To establish that this second problem can theoretically have a solution, consider solving (A11) for two modified pressure gradients
$\boldsymbol{\tilde {G}}_A$
and
$\boldsymbol{\tilde {G}}_B$
. Let
$\boldsymbol{\tilde {u}}_A$
and
$\boldsymbol{\tilde {u}}_B$
be the corresponding solutions. We can see that
$\boldsymbol{\tilde {u}}_A$
is a test function for
$\boldsymbol{\tilde {u}}_B$
and
$\boldsymbol{\tilde {u}}_B$
is a test function for
$\boldsymbol{\tilde {u}}_A$
. On substituting into (A11) for both the
$\boldsymbol{\tilde {u}}_A$
and
$\boldsymbol{\tilde {u}}_B$
problems, summing the inequalities and rearranging, we find that
\begin{eqnarray} [\boldsymbol{\tilde {G}}_B - \boldsymbol{\tilde {G}}_A] \boldsymbol{\cdot }[ \langle \boldsymbol{\tilde {u}}_B\rangle -\langle \boldsymbol{\tilde {u}}_A \rangle ] & \geqslant & a(\boldsymbol{\tilde {u}}_A, \boldsymbol{\tilde {u}}_A-\boldsymbol{\tilde {u}}_B) - a(\boldsymbol{\tilde {u}}_B, \boldsymbol{\tilde {u}}_A-\boldsymbol{\tilde {u}}_B) \geqslant 0 . \end{eqnarray}
The last inequality in (A13) follows from the monotonicity of the operator
$a$
that is the sub-differential of a convex functional. The inequality is strict unless
$\boldsymbol{\tilde {u}}_A = \boldsymbol{\tilde {u}}_B = 0$
, which can happen for different
$\boldsymbol{\tilde {G}}_A$
and
$\boldsymbol{\tilde {G}}_B$
, in the case that the yield stresses are large enough to prevent motion. The positivity of
$[\boldsymbol{\tilde {G}}_B - \boldsymbol{\tilde {G}}_A] \boldsymbol{\cdot }[ \langle \boldsymbol{\tilde {u}}_B\rangle -\langle \boldsymbol{\tilde {u}}_A \rangle ]$
suggests that one can simply increase
$\boldsymbol{\tilde {G}}$
in the right direction to either increase/decrease
$\langle \boldsymbol{\tilde {u}} \rangle$
towards a target mean velocity.
A.2. Augmented Lagrangian method
A common solution method for the minimisation problem (A12) involves a relaxation of
${{\rm d} \boldsymbol{\tilde {v}}}/{{\rm d} \tilde {y}} \longmapsto \boldsymbol{q}$
in the definition of
$J(\boldsymbol{\tilde {v}})$
, which removes the non-differentiability of
$j(\boldsymbol{\tilde {v}})$
. The constraint
${{\rm d} \boldsymbol{\tilde {v}}}/{{\rm d} \tilde {y}} = \boldsymbol{q}$
is then enforced by a Lagrange multiplier and the minimisation is usually also stabilised by augmenting with an additional penalty term. This results in the following saddle point problem:
\begin{equation} \max _{\boldsymbol{\lambda }}, \min _{\boldsymbol{\tilde {v}},\boldsymbol{q}} \mathcal{L}_r(\boldsymbol{\tilde {v}},\boldsymbol{q},\boldsymbol{\lambda }). \end{equation}
Here the functional
$\mathcal{L}_r$
is as follows:
\begin{align} \mathcal{L}_r(\boldsymbol{\tilde {v}},\boldsymbol{q},\boldsymbol{\lambda }) & = \frac { \tilde {\kappa }_1 \left\langle |\boldsymbol{q}|^{n_1+1} \right\rangle _1 }{n_1+1} + \frac { \tilde {\kappa }_2 \left\langle |\boldsymbol{q}|^{n_2+1} \right\rangle _2 }{n_2+1} + \tilde {\unicode{x03C4}}_{1,Y} \langle |\boldsymbol{q}| \rangle _1 + \tilde {\unicode{x03C4}}_{2,Y} \langle |\boldsymbol{q}| \rangle _2 - \boldsymbol{\tilde {G}} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle \nonumber \\ &\quad + (1-\bar {c}) \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle _2 - \bar {c} \boldsymbol{\tilde {G}_b} \boldsymbol{\cdot }\langle \boldsymbol{\tilde {v}} \rangle _1 + \left\langle \boldsymbol{\lambda } \boldsymbol{\cdot }\left( \frac {{\rm d} \boldsymbol{\tilde {v}}}{{\rm d} \tilde {y}} - \boldsymbol{q} \right) \right\rangle \nonumber \\ &\quad + \frac {r}{2}\left\langle \left(\frac {{\rm d}\boldsymbol{\tilde {v}}}{{\rm d}y}-\boldsymbol{q}\right) \boldsymbol{\cdot }\left(\frac {{\rm d}\boldsymbol{\tilde {v}}}{{\rm d}y}-\boldsymbol{q}\right) \right\rangle . \end{align}
This is resolved iteratively using an Uzawa algorithm, in which new values of
$\boldsymbol{\tilde {u}}$
,
$\boldsymbol{q}$
and
$\boldsymbol{\lambda }$
are found, independently at each iteration.
A.2.1. Fixed pressure gradient problem
We first treat the problem in which
$\boldsymbol{\tilde {G}}$
and
$\boldsymbol{\tilde {G}_b}$
are specified. At the
$(k+1)$
th iteration,assuming that
$\boldsymbol{q}^{k}$
and
$\boldsymbol{\lambda }^{k}$
are known from the previous iteration, the optimality condition for
$\boldsymbol{\tilde {v}}$
, which defines
$\boldsymbol{\tilde {u}}^{k+1}$
, is given by
\begin{align} r\frac {d^2\boldsymbol{\tilde {u}}^{k+1}}{{\rm d}\tilde {y}^2} &= r\frac {{\rm d}\boldsymbol{q}^{k}}{{\rm d}\tilde {y}}-\frac {{\rm d}\boldsymbol{\lambda }^{k}}{{\rm d}\tilde {y} }-\boldsymbol{\tilde {G}} + (1-\bar {c}) \boldsymbol{\tilde {G}_b}, \quad \tilde {y} \in [0, \bar {c}),\nonumber\\ r\frac {d^2 \boldsymbol{\tilde {u}}^{k+1}}{{\rm d}\tilde {y}^2} &= r\frac {{\rm d}\boldsymbol{q}^{k}}{{\rm d}\tilde {y}}-\frac {{\rm d}\boldsymbol{\lambda }^{k}}{{\rm d}\tilde {y}}-\boldsymbol{\tilde {G}}-\bar {c}\boldsymbol{\tilde {G}_b}, \quad \tilde {y} \in (\bar {c}, 1], \end{align}
with conditions of symmetry at
$\tilde {y} =0$
, continuity at
$\tilde {y} = \bar {c}$
and
$\boldsymbol{\tilde {u}}^{k+1}=0$
at
$\tilde {y} =1$
.
Since we expect that the sequence
$\boldsymbol{\lambda }^k$
converges to
$\boldsymbol{\tilde {\unicode{x03C4}}}$
, which varies linearly with
$\tilde {y}$
, we split
$\boldsymbol{\lambda }^k = \boldsymbol{\lambda }_0 + \boldsymbol{\tilde {\lambda }}^k$
, where we select
$\boldsymbol{\lambda }_0$
to be piecewise linear and balance the
$\boldsymbol{\tilde {G}}$
and
$\boldsymbol{\tilde {G}_b}$
terms:
\begin{equation} \boldsymbol{\lambda }_0(\tilde {y}) = \left \{ \!\!\begin{array}{lcl} - \boldsymbol{\tilde {G}}\tilde {y} + (1-\bar {c}){\boldsymbol{\tilde {G}_b}}\tilde {y}, &\,\,\,& \tilde {y} \in [0, \bar {c}), \\ - \boldsymbol{\tilde {G}}\tilde {y} + (1-\tilde {y}){\boldsymbol{\tilde {G}_b}} \bar {c}, & & \tilde {y} \in (\bar {c},1]. \end{array} \right . \end{equation}
Using (A17), we now simplify (A16) to
\begin{equation} r\frac {d^2\boldsymbol{\tilde {u}}^{k+1}}{{\rm d}\tilde {y}^2} = r\frac {{\rm d}\boldsymbol{q}^{k}}{{\rm d}\tilde {y}}-\frac {{\rm d} \boldsymbol{\tilde {\lambda }}^{k}}{{\rm d}\tilde {y}},\,\,\,\, \Rightarrow \,\,\,\, r\frac {{\rm d}\boldsymbol{\tilde {u}}^{k+1}}{{\rm d}\tilde {y}} = r \boldsymbol{q}^{k}(\tilde {y}) - \boldsymbol{\tilde {\lambda }}^{k}(\tilde {y}) - r \boldsymbol{q}^{k}(0) + \boldsymbol{\tilde {\lambda }}^{k}(0) . \end{equation}
This is integrated from the wall, where
$\boldsymbol{\tilde {u}}^{k+1}(1) = 0$
.
To find
$\boldsymbol{q}^{k+1}$
, we may minimise locally, i.e.
\begin{equation} \boldsymbol{K}(\boldsymbol{q}) = \frac { \tilde {\kappa }_j |\boldsymbol{q}|^{n_j+1} }{n_j+1} + \frac { r |\boldsymbol{q}|^{2} }{2} + \tilde {\unicode{x03C4}}_{Y,j} |\boldsymbol{q}| - \boldsymbol{m} \boldsymbol{\cdot }\boldsymbol{q}, \end{equation}
where
$\boldsymbol{m} = \boldsymbol{\lambda } + \boldsymbol{\tilde {\lambda }}^k + r ({{\rm d}\boldsymbol{\tilde {u}}^{k+1}}/{{\rm d}y})$
and
$j = 1,2$
according to which fluid layer we solve in.
The optimality conditions for
$\boldsymbol{q}^{k+1}$
are
\begin{equation} \boldsymbol{q}^{k+1} = \left \{ \!\!\begin{array}{lcl} 0, & \,\,\,\, & |\boldsymbol{m}|\leqslant \unicode{x03C4}_{Y,j}, \\ \theta \boldsymbol{m}, & & |\boldsymbol{m}| \gt \unicode{x03C4}_{Y,j}. \\ & & \mbox{where } \tilde {\kappa }_j \theta ^{n_j} |\boldsymbol{m}|^{n_j} + r \theta |\boldsymbol{m}| = |\boldsymbol{m}|-\unicode{x03C4}_{Y,j} , \end{array} \right . \end{equation}
which requires iterative solution of an equation of one variable to find
$\theta$
. Finally, we update
$\tilde {\boldsymbol{\lambda }}$
:
\begin{equation} \boldsymbol{\tilde {\lambda }}^{k+1}=\boldsymbol{\tilde {\lambda }}^{k} + \rho (\frac {{\rm d}\boldsymbol{u}^{k+1}}{{\rm d}y}-\boldsymbol{q}^{k+1}) . \end{equation}
The positive constant parameters
$r$
and
$\rho$
represent the numerical parameters of the Uzawa algorithm and can be selected to guarantee convergence (Fortin & Glowinski Reference Fortin and Glowinski2000; Glowinski Reference Glowinski2013).
A.2.2. Fixed average velocity
In the second closure model scenario it is necessary to calculate the pressure gradient
$\boldsymbol{\tilde {G}}$
for an imposed average velocity
$\boldsymbol{\bar {u}}^*$
and buoyancy
$\boldsymbol{\tilde {G}_b}$
. There are a number of ways in which this may be solved. For example, one could guess
$\boldsymbol{\tilde {G}}^k$
and then use the monotonicity result (A13) to update
$\boldsymbol{\tilde {G}}^{k+1}$
in the direction of the discrepancy between
$\boldsymbol{\bar {u}}^*$
and the computed
$\boldsymbol{\bar {u}}^k$
.
Instead, our underlying approach is to use the saddle point problem (A14), but exploit the linearity of the first step to satisfy the flow rate constraint by assigning a new value
$\boldsymbol{\tilde {G}}^{k+1}$
tied to the velocity solution at each iteration. To start, let us define
$\boldsymbol{\bar {u}}_{P,0}$
as the solution of the following plane channel Newtonian Poiseuille flow:
\begin{equation} r\frac {d^2\boldsymbol{u}_{P,0}}{{\rm d}y^2} = -\boldsymbol{\tilde {G}}_0, \qquad\frac {{\rm d}\boldsymbol{u}_{P,0}}{{\rm d}y}(0) = 0,\qquad \boldsymbol{u}_{P,0}(1) = 0. \end{equation}
Here
$\boldsymbol{\tilde {G}}_0 = (1,1)$
, i.e. unit pressure gradient. This may be solved analytically, but generally we solve computationally using the same discretisation as for the velocities below. Having solved, we also calculate the average velocity
$\bar {\boldsymbol{u}}_{P,0}$
, which analytically is given by
$\bar {\boldsymbol{u}}_{P,0} = (1/3r,1/3r)$
. We also define
$\boldsymbol{\lambda }_0$
to balance the buoyancy
$\boldsymbol{\tilde {G}_b}$
:
\begin{equation} \boldsymbol{\lambda }_0(\tilde {y}) = \left \{ \!\!\begin{array}{lcl} (1-\bar {c}){\boldsymbol{\tilde {G}_b}}\tilde {y}, &\,\,\,& \tilde {y} \in [0, \bar {c}), \\ (1-\tilde {y}){\boldsymbol{\tilde {G}_b}} \bar {c}, & & \tilde {y} \in (\bar {c},1]. \end{array} \right . \end{equation}
At the
$(k+1)$
th iteration for (A14), assuming that
$\boldsymbol{q}^{k}$
and
$\boldsymbol{\lambda }^{k}$
are known from the previous iteration, the optimality condition for
$\boldsymbol{\tilde {v}}$
is given by
\begin{align} r\frac {d^2\boldsymbol{\tilde {u}}^{k+1}}{{\rm d}\tilde {y}^2} &= r\frac {{\rm d}\boldsymbol{q}^{k}}{{\rm d}\tilde {y}}-\frac {{\rm d}\boldsymbol{\lambda }^{k}}{{\rm d}\tilde {y} }-\boldsymbol{\tilde {G}}^{k+1} + (1-\bar {c}) \boldsymbol{\tilde {G}_b}, \quad \tilde {y} \in [0, \bar {c}),\nonumber\\ r\frac {d^2 \boldsymbol{\tilde {u}}^{k+1}}{{\rm d}\tilde {y}^2} &= r\frac {{\rm d}\boldsymbol{q}^{k}}{{\rm d}\tilde {y}}-\frac {{\rm d}\boldsymbol{\lambda }^{k}}{{\rm d}\tilde {y}}-\boldsymbol{\tilde {G}}^{k+1}-\bar {c}\boldsymbol{\tilde {G}_b}, \quad \tilde {y} \in (\bar {c}, 1], \end{align}
again with symmetry at
$\tilde {y} =0$
, continuity conditions at
$\tilde {y} = \bar {c}$
and
$\boldsymbol{\tilde {u}}^{k+1}=0$
at
$\tilde {y} =1$
. Now, however, we must solve this such that
$\boldsymbol{\bar {\tilde {u}}}^{k+1} = \boldsymbol{\bar {u}}^*$
. To do this we split
$\boldsymbol{\tilde {u}}^{k+1}$
, i.e.
\begin{align} \boldsymbol{\tilde {u}}^{k+1} = \boldsymbol{\tilde {u}}_I^{k+1} + \boldsymbol{u}_P^{k+1}, \end{align}
where
$\boldsymbol{\tilde {u}}_I^{k+1}$
takes care of the inhomogeneities in (A24) and
$\boldsymbol{u}_P^{k+1}$
ensures that the mean velocity condition is satisfied.
More concretely, we again define
$\boldsymbol{\tilde {\lambda }}^k$
as
$\boldsymbol{\tilde {\lambda }}^k = \boldsymbol{\lambda }^k - \boldsymbol{\lambda }_0$
and find
$\boldsymbol{\tilde {u}}_I$
from
\begin{equation} r\frac {d^2\boldsymbol{\tilde {u}}_I^{k+1}}{{\rm d}\tilde {y}^2} = r\frac {{\rm d}\boldsymbol{q}^{k}}{{\rm d}\tilde {y}}-\frac {{\rm d} \boldsymbol{\tilde {\lambda }}^{k}}{{\rm d}\tilde {y}}\,\,\,\, \Rightarrow \,\,\,\, r\frac {{\rm d}\boldsymbol{\tilde {u}}_I^{k+1}}{{\rm d}\tilde {y}} = r \boldsymbol{q}^{k}(\tilde {y}) - \boldsymbol{\tilde {\lambda }}^{k}(\tilde {y}) - r \boldsymbol{q}^{k}(0) + \boldsymbol{\tilde {\lambda }}^{k}(0) . \end{equation}
This is integrated from the wall, where
$\boldsymbol{\tilde {u}}_I^{k+1}(1) = 0$
. We now average over
$[0,1]$
to give
$\boldsymbol{\bar {\tilde {u}}}_I^{k+1}$
. Therefore, we must select
$\boldsymbol{u}_P^{k+1}$
in order that
\begin{equation} \bar {\boldsymbol{u}}_P^{k+1} = \boldsymbol{\bar {u}}^* - \boldsymbol{\bar {\tilde {u}}}_I^{k+1}. \end{equation}
This is accomplished by taking
$\bar {\boldsymbol{u}}_P^{k+1} = (\alpha _\phi u_{P,0,\phi },\alpha _\xi u_{P,0,\xi })$
, where
\begin{align} \alpha _\phi = \frac {\bar {u}^*_{\phi } - \bar {\tilde {u}}_{I,\phi }^{k+1}}{\bar {u}_{P,0,\phi }}, \qquad\alpha _\xi = \frac {\bar {u}^*_{\xi } - \bar {\tilde {u}}_{I,\xi }^{k+1}}{\bar {u}_{P,0,\xi }} , \end{align}
and
$\boldsymbol{\tilde {G}}^{k+1} = (\alpha _\phi ,\alpha _\xi )$
.
The optimality conditions for
$\boldsymbol{q}^{k+1}$
are again
\begin{equation} \boldsymbol{q}^{k+1} = \left \{ \!\!\begin{array}{lcl} 0, & \,\,\,\, & |\boldsymbol{m}|\leqslant \unicode{x03C4}_{Y,j}, \\[6pt] \theta \boldsymbol{m}, & & |\boldsymbol{m}| \gt \unicode{x03C4}_{Y,j}. \\[6pt] & & \mbox{where } \tilde {\kappa }_j \theta ^{n_j} |\boldsymbol{m}|^{n_j} + r \theta |\boldsymbol{m}| = |\boldsymbol{m}|-\unicode{x03C4}_{Y,j} , \end{array} \right . \end{equation}
where
$\boldsymbol{m} = \boldsymbol{\lambda } + \boldsymbol{\tilde {\lambda }}^k + r ({{\rm d}\boldsymbol{\tilde {u}}^{k+1}}/{{\rm d}y})$
and
$j = 1,2$
according to which fluid layer we solve in. Note that the pressure gradient
$\boldsymbol{\tilde {G}}^{k+1}$
now comes into this via
$\bar {\boldsymbol{u}}_P^{k+1}$
whereas buoyancy is represented via
$\boldsymbol{\lambda }_0$
. Finally, we update
$\tilde {\boldsymbol{\lambda }}$
:
\begin{equation} \boldsymbol{\tilde {\lambda }}^{k+1}=\boldsymbol{\tilde {\lambda }}^{k} + \rho \left(\frac {{\rm d}\boldsymbol{\tilde {u}}^{k+1}}{{\rm d}y}-\boldsymbol{q}^{k+1}\right) . \end{equation}
The overall iteration is thus quite similar to the fixed pressure gradient problem.
The norms of
$\boldsymbol{\tilde {u}}, \boldsymbol{\tilde {q}},\text{ and } \boldsymbol{\tilde {\lambda }}$
are plotted against the number of iterations (
$k$
) for Newtonian, shear-thinning and Herschel–Bulkley fluids for isodense flows (
$b=0$
) at
$\kappa _1=\kappa _2=1$
. Here,
$\bar {c}=0.5$
and the mesh size is set to
$N_Y=100$
. The broken blue line plotted in (a) represents the stopping tolerance for these computations.
A.2.3. Discretisation and illustrations of convergence
We discretise the equations using a finite difference discretisation with a staggered regular rectangular mesh, i.e. we divide the width of the channel into approximately
$N_y$
mesh points with
$\approx N_y(1-\bar {c})$
for the displaced fluid and
$\approx N_y\bar {c}$
for the displacing fluid. We define the discretised function
$\boldsymbol{\tilde {u}}^k$
at the cell edges and
$\boldsymbol{\lambda }^k, \boldsymbol{q}^k$
,
$\bar {c}$
at the cell centre. This definition is useful for calculating the gradients arising in the equations (for both fixed
$\boldsymbol{G}$
and
$\boldsymbol{\bar {u}}$
closures).
We use a simple bisection method for calculating
$\boldsymbol{q}^{k+1}$
and the constant numerical parameters of the Uzawa algorithm (
$\rho$
,
$r$
) are fixed at
$\rho =1$
,
$r=1$
that satisfies the convergence condition:
\begin{equation} \rho \lt \frac {1+\sqrt {5}}{2}r \end{equation}
(Fortin & Glowinski Reference Fortin and Glowinski2000). We then monitor the norms,
$||\boldsymbol{\tilde {u}}^{k+1}-\boldsymbol{\tilde {u}}^{k}||_p$
,
$||\boldsymbol{\tilde {q}}^{k+1}-\boldsymbol{\tilde {q}}^{k}||_p$
and
$||\boldsymbol{\tilde {\lambda }}^{k}||_p$
, iterating until a desired tolerance value is achieved. Here
$||.||_p$
corresponds to the
$p$
norm of the variables, with
$p=1+\min \{n_1,n_2\}$
.
Figures 18 and 19 illustrate the convergence of the algorithm. For two identical fluids, from figure 18(a–c), it can be seen that increasing the power-law index (
$n$
) improves the convergence as does reducing the yield stress, i.e. the Newtonian flow converges fastest.
In figure 19 the impact of the buoyancy number (
$b$
) and the viscosity ratio (
$m$
) on the convergence is explored. For the isodense flows, the convergence rate increases when the displaced fluid is more viscous than the displacing fluid; see figure 19(a). The convergence also improves when the density of the displacing fluid becomes higher than the displaced fluid, as shown in figure 19(b). Comparing figure 19(c) with figure 19(b), the negative
$b$
slows convergence, where buoyancy has a destabilising effect on the flow.
The norm of
$\boldsymbol{\tilde {u}}$
is plotted against the number of iterations (
$k$
) for Newtonian, shear-thinning and Herschel–Bulkley fluids at (a)
$b=0$
(isodensity), (b)
$b=15$
(displacing fluid heavier than displaced) and (c)
$b=-15$
(displaced fluid heavier than displacing). Here,
$\bar {c}=0.5$
and the mesh size is set to
$N_Y=100$
. The black lines represent the tolerance of this computation.
The other physical parameter that has an impact on the convergence rate is the initial guess for
$\boldsymbol{\tilde {G}}$
and
$\boldsymbol{\tilde {q}}$
. There can also be sensitivity to the numerical parameters (
$\rho$
,
$r$
), which we have not explored.
In choosing a tolerance, our integration method is second order in the mesh size
$1/N_y$
, so tolerances much smaller than
$1/N_y^2$
do not improve the physical representation. Even though the iterations here are relatively fast, if one is considering using this solution method to calculate the D2DGA closure expressions as part of the (no doubt iterative) solution of an evolution problem for the streamfunction and concentration(s) in a long wellbore domain, it is not feasible. Instead, it is advisable to pre-compute the closure expressions needed for a full range of reduced rheological parameters and over restricted ranges of
$b$
,
$\boldsymbol{\tilde {G}}$
, etc.
A.3. Examples of axial velocity profiles
Figure 20 illustrates various axial velocity profiles (
$w(y)$
) obtained using the fixed average velocity closure, i.e. the pressure gradient is found with the student. Note that in figure 20(a–e) the modified pressure gradient is higher than the buoyancy number (
$G_\xi \gt b$
), and accordingly, fluids are pushed upwards (in the positive direction). Flow reversal is displayed in figure 20( f) under the condition of
$b\lt 0$
, where the displaced fluid is denser than the displacing fluid.
Newtonian fluid flows at
$\kappa _1=\kappa _2=1$
(a), viscoplastic fluid flows with
$n_1=n_2=0.5$
at (b)
$\tau _{Y,1}=0,\,\tau _{Y,2}=0.89,\,\kappa _1=\kappa _2=0.105$
, (c)
$\tau _{Y,1}=0.123,\,\tau _{Y,2}=0.98,\,\kappa _1=\kappa _2=0.014$
, (d)
$\tau _{Y,1}=0.77,\,\tau _{Y,2}=0.038,\,\kappa _1=\kappa _2=0.227$
, (e)
$\tau _{Y,1}=0.98,\,\tau _{Y,2}=0.123,\,\kappa _1=\kappa _2=0.014$
, ( f)
$\tau _{Y,1}=0,\,\tau _{Y,2}=0.89,\,\kappa _1=\kappa _2=0.105$
.
Open access
