2. Methodology

As discussed in § 1 the main objectives of this work are to expose deficiencies in 2-D gap-averaged displacement models of primary cementing and to seek a viable way to improve them. As always in practical fluid mechanics, the choice of model is a balance between accuracy and computational cost. Here we introduce the different models used in our study.

2.1. The 2DGA model

We outline here the current 2DGA model for annular displacement flows, with more detail later in § 4.1, where we improve the model. We adopt the formulation of Maleki & Frigaard (Reference Maleki and Frigaard2017), where the model is more comprehensive than that earlier (Bittleston et al. Reference Bittleston, Ferguson and Frigaard2002). The geometry of the narrow eccentric annulus is mapped and unwrapped into a Hele-Shaw cell, see figure 2. The key assumption of the model is that the mean annular gap ($2\hat {d}$) is much narrower compared with both the mean circumference ($2 {\rm \pi}\hat {r}_a^*$), and to the length scale of any axial geometric variations. The narrow gap approximation is that $\hat {d}/({\rm \pi} \hat {r}_a^*) = \delta /{\rm \pi} \ll 1$. The velocity and fluid concentration profiles are then averaged in the radial direction, which reduces the problem to a 2-D one. The length scale assumptions also allow us to assume the flow is locally a shear flow. With the additional assumption that the gap-averaged quantities describe the fluid present, we are able to solve the 1-D shear flow in the direction of the modified pressure gradient locally. This provides the closure expressions that are employed in the 2DGA model. The model consists of a nonlinear elliptic equation for the stream function and a sequence of advection equations for the fluid concentrations.

Figure 2. Geometry of the narrow eccentric annulus mapped to the Hele-Shaw cell, for the 2DGA model.

All dimensional variables are with the ‘hat’ accent, e.g. $\hat {w}_0$, and dimensionless variables without. The governing dimensionless equations are listed as follows in (2.1)–(2.3), we discuss only laminar flow regimes in this paper:

(2.1)\begin{gather} \frac{\partial}{\partial t}[Hr_{a}\bar{c}_{k}]+\frac{\partial}{\partial \phi}[H\bar{v}\bar{c}_{k}]+\frac{\partial}{\partial \xi}[r_a H\bar{w}\bar{c}_{k}]=0,\quad k=1,2,\ldots,K, \end{gather}

(2.2)\begin{gather} \boldsymbol{\nabla}_{a} \boldsymbol{\cdot} ( 2r_a H\bar{v},2r_a H\bar{w}) = 0, \quad \Rightarrow\quad \frac{\partial \varPsi }{\partial \phi }=2Hr_{a}\bar{w}, \quad \frac{\partial \varPsi }{\partial \xi }={-}2H\bar{v}, \end{gather}

(2.3)\begin{gather} {\boldsymbol{\nabla}_{a}\boldsymbol{\cdot} [ \boldsymbol{S} + \boldsymbol{b}] = 0.} \end{gather}

Equation (2.1) is a transport equation for the fluid concentrations $\bar {c}_{k}$ for $k=1,2,\ldots,K$, which are advected with the gap-averaged velocity $(\bar {v},\bar {w})$. The scaled half-gap-width is $H(\phi,\xi )$ and the mean radius at depth $\xi$ is $r_a(\xi )$. In laminar miscible flows the effects of molecular diffusion are minimal on the time scales and length scales of relevance. Thus, diffusive effects are absent on the right-hand side of (2.1), but not for turbulent regimes, where Taylor dispersion dominates (Maleki & Frigaard Reference Maleki and Frigaard2017). In deriving (2.1), azimuthal and axial lengths have been scaled with a representative (half-)circumference of the annulus ${\rm \pi} \hat {r}_a^*$. Velocities have been scaled with a mean pumping velocity $\hat {w}_0$, and times with ${\rm \pi} \hat {r}_a^*/\hat {w}_0$. The mean radius $\hat {r}_a^*$ is defined by averaging along the annular flow path. The gap averaged velocity field is incompressible, allowing it to be expressed in terms of a stream function $\varPsi$, i.e. (2.2). The operator ${\nabla }_{a}$ is a radial divergence operator, i.e.

(2.4)\begin{equation} {\nabla}_{a} = \left( \frac{1}{r_{a}}\frac{\partial }{\partial \phi }, \frac{\partial }{\partial \xi } \right), \end{equation}

with ${\nabla }_{a}$ is the corresponding gradient operator.

The stream function satisfies the elliptic equation (2.3), where

(2.5a,b)\begin{equation} \boldsymbol{S} = \left[\frac{r_{a}\tau_w (\left | {\nabla}_{a} \varPsi \right |)}{H \left | {\nabla}_{a}\varPsi \right |} \right] {\nabla}_{a}\varPsi ,\quad \boldsymbol{b} = \frac{r_a (\rho-1)}{{F^2}} (\cos \beta, \sin {\rm \pi}\phi \sin \beta ) . \end{equation}

The function $\boldsymbol {S}$ represents a modified pressure gradient, which we see has magnitude $r_a \tau _w ( | {\nabla }_{a} \varPsi |)/H$. Here $\tau _w$ is the wall shear stress. These quantities have been made dimensionless using a stress scale $\hat {\tau }_0$. In the definition of $\boldsymbol {b}$, $\rho$ is the density (scaled with the density of the first fluid), and $F$ is the Froude number,

(2.6)\begin{equation} {F} = \sqrt{\frac{\hat{\tau}_{0}}{ \hat{\rho}_1 \hat{g} \delta_0 \hat{r}_{a}^* }} . \end{equation}

For clarity, here we fix $\hat {\tau }_{0}=\hat {\mu }_1 \hat {w}_0/\hat {d}$, i.e. based on a viscous stress in the displaced fluid. However, more generally $\hat {\tau }_{0}$ could be chosen to represent any relevant stress (Maleki & Frigaard Reference Maleki and Frigaard2017), depending on the dominant flow regime and rheology. Thus, we see that $\boldsymbol {b}$ represents the effects of buoyancy in this model; together with the imposed flow, buoyancy gradients strongly influence the flow.

The dependency of $\tau _w$ on the local flow rate (mean velocity), on the fluid concentrations and on the rheological parameters will depend on the specific fluid properties, i.e. the model is completed with a closure expression for $\tau _w$. Generally speaking, in laminar flow regimes cementing involves shear thinning, yield-stress fluids. The yield stress of these fluids $\tau _Y$ plays an important role industrially, defining a limiting pressure gradient $r_{a}\tau _{Y}/H$, that must be exceeded locally in order for the fluids to flow; otherwise $| {\nabla }_{a} \varPsi |=0$ and the fluids are stuck. The 2-D model was first analysed by Pelipenko & Frigaard (Reference Pelipenko and Frigaard2004a,Reference Pelipenko and Frigaardb) and there are a large number of works that use this underlying framework.

Numerically, we follow the procedures described in Maleki & Frigaard (Reference Maleki and Frigaard2017). The equations are discretised on a regular rectangular mesh, with staggered mesh for the stream function versus concentration. We see that the flow evolves in time only via (2.1), which is time-advanced explicitly using a flux corrected transport scheme. At each time step the nonlinear elliptic stream function equation (2.3), is solved using an iterative augmented Lagrangian algorithm. The velocities are constructed from $\varPsi$ via (2.2) and used for the fluid concentrations in (2.1).

2.2. The 3-D model

The 3-D numerical simulations are computed using a finite volume method within the open source computational fluid dynamics toolbox OpenFOAM (http://www.openfoam.com). The code used here is an updated version of that used in Etrati & Frigaard (Reference Etrati and Frigaard2018) for studying miscible displacement flows in pipe flows, wherein the results were successfully compared with experiments. More recently, the code has been used by Sarmadi et al. (Reference Sarmadi, Renteria and Frigaard2021), comparing with the experimental results of Renteria & Frigaard (Reference Renteria and Frigaard2020). Briefly, the full Navier–Stokes equations are solved, coupled to an advection equation for the concentration of fluid 2, which is the displacing fluid. The volume of fluid method is used to capture the interface between the fluids. The momentum, continuity and concentration equations are

(2.7)\begin{gather} \hat{\rho}\left(\frac{\partial \hat{\boldsymbol{u}}}{\partial \hat{t}}+\hat{\boldsymbol{u}}\boldsymbol{\cdot} \hat{\nabla} \hat{\boldsymbol{u}} \right)={-}\hat{\nabla} \hat{p}+\hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\hat{\boldsymbol{\tau}}+\hat{\rho}\hat{\boldsymbol{g}}, \end{gather}

(2.8)\begin{gather} \hat{\boldsymbol{\nabla}}\boldsymbol{\cdot}\hat{\boldsymbol{u}}=0, \end{gather}

(2.9)\begin{gather} \frac{\partial c}{\partial \hat{t}}+\hat{\boldsymbol{u}} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} c=0. \end{gather}

Here, $\hat {\boldsymbol {u}}$ and $\hat {p}$ are the velocity and pressure; $\hat {\boldsymbol {\tau } }=\hat {\mu }\hat {\dot {\boldsymbol {\gamma } }}(\hat {\boldsymbol {u}})$ is the deviatoric stress tensor, with $\hat {\dot {\boldsymbol {\gamma } }}$ the rate of strain tensor. In this paper we only consider two Newtonian fluids. The concentration $c\in [0,1]$ enters via the fluid properties: $\hat {\rho }$ and $\hat {\mu }$ are the mixture density and viscosity, defined by

(2.10)\begin{gather} \hat{\rho} =(1-c)\hat{\rho}_{1}+c\hat{\rho}_{2}, \end{gather}

(2.11)\begin{gather} \hat{\mu}=\hat{\mu} _{1}\exp\left[c\,\mathrm{log}\left(\frac{\hat{\mu} _{2}}{\hat{\mu} _{1}}\right)\right]. \end{gather}

The subscript ‘1’ represents the displaced fluid and ‘2’ represents the displacing fluid. We have also used linear interpolation for the viscosity closure, with only mild differences.

A uniform velocity is imposed at the inlet, no-slip is applied at the annulus walls and an outflow condition set at the outlet. The pressure is fixed at the outlet. The PIMPLE algorithm is used to solve the momentum equation, and the implicit second-order Crank–Nicolson method is employed for the time discretisation. Mesh sensitivity/convergence results are given in Sarmadi et al. (Reference Sarmadi, Renteria and Frigaard2021). Figure 3 illustrates the computational mesh used for one of our cases with eccentricity $e = 0.4$, with $20 \times 100 \times 400$ cells in the radial, azimuthal and axial directions. Figure 3(a) shows half of the annulus, taken from the midplane ($x = 0$), and figure 3(b) shows a middle slice surface. The mesh is refined close to the wall to capture residual wall layers, as common in displacement flows. The computations are carried out in parallel on a multiprocessor machine, generally using 24 cores.

Figure 3. Mesh resolution in (a) axial (vertical) direction and (b) the cross-section for an eccentric case with $e = 0.4$.

2.3. Scope of the study

The aim is to study and improve the 2DGA model, making use of the 3-D computations. We also later wish to use these same models for the purposes of experimental comparison. Thus, we set the annulus dimensions and other parameters of our study based on the experimental set-up of Renteria & Frigaard (Reference Renteria and Frigaard2020), which we have also used for our experimental study. The total length of the annulus is 4.8 m, the outer and inner radii are 22.23 mm and 17.46 mm, respectively. This gives a scaled length $Z \approx 62$ for the 2DGA model, i.e. the length is $\approx 62 \times$ the mean circumferential length scale ${\rm \pi} \hat {r}_a^*$, (and is $>1000 \times$ the mean annular gap width). The eccentricity can be adjusted from 0 to 1. Unlike Renteria & Frigaard (Reference Renteria and Frigaard2020) and Sarmadi et al. (Reference Sarmadi, Renteria and Frigaard2021), here we study vertical annuli.

For two Newtonian fluids the flow depends on at least seven dimensionless parameters, and more if the fluids are non-Newtonian. Having adopted $\hat {\tau }_{0}=\hat {\mu }_1 \hat {w}_0/\hat {d}$, as our (viscous) stress scale, the key dimensionless parameters can be taken as the Reynolds number, $Re=\hat {\rho }_{1}\hat {w}_{0}\hat {d}/\hat {\mu }_{1}$, the viscosity ratio, $m=\hat {\mu } _{1}/\hat {\mu } _{2}$, the buoyancy number, $b=(\hat {\rho }_{2}-\hat {\rho }_{1} )\hat {g}\hat {d}^{2}/(\hat {\mu }_{1}\hat {w}_{0})$, a density ratio (or Atwood number), the pipe inclination and eccentricity. We see that $b$ represents the size of the buoyancy vector $\boldsymbol {b}$ in the 2DGA model. Rather than explore all parameters systematically, we instead study a selection of 10 cases that expose key features of vertical displacement. The cases and corresponding parameters are listed in tables 1 and 2. For simplicity, all fluid pairs are Newtonian and the annulus is vertical. Each simulation case is run for a time equivalent to two annular volumes of pumped fluid, using both styles of model.

Table 1. Fluid properties for the designed simulation cases.

Table 2. Corresponding dimensionless parameters for the designed simulation cases.

2.4. Mass conservation properties

To verify that both models are correctly programmed and solved, we select cases 2, 4, 6 and 8, each of which has different flow rates. We integrate the equations forward and compute the total volume of the displacing fluid in the annulus $(\hat {V} _{2})$ at successive time steps. This is plotted in figure 4, in which we compare the two computations with the theoretical volume (from the inflow rate). Stopping the simulations at any particular time and dividing by the annulus volume, gives a fraction of the annulus displaced, i.e. the displacement efficiency. We see that in each case, the 2DGA and 3-D models compute near-identical volumes up until the breakthrough time, i.e. that time at which the displacing fluid first exits the annulus. The breakthrough time is generally earlier for the 3-D simulation, due to dispersion ahead of the main front. We also observe earlier breakthrough for cases 6 and 8, for both of which the displacing fluid is less viscous than the displaced fluid ($m > 1$). The same verification is made for the other cases, with similar agreement; not shown for reasons of brevity.

Figure 4. Verification of volume conservation for the 2DGA and 3-D models: (a) cases 2 and 4; (b) cases 6 and 8. The black lines represent the total annulus volume. The solid lines with other colours represent the theoretical value based on the constant flow rate. The 2DGA and 3-D results are marked by different symbols, as indicated.

3. Successes and failures of the Hele-Shaw approach

Now we explore some of these cases in more detail, here and in the Appendix. We start with case 4, which is representative of a successful displacement flow, (commonly found if $e$ is modest and $b$ is large). Figure 5 shows the comparison of 2-D gap-averaged concentration of case 4 at successive times as indicated. The variable $\phi \in [0,1]$ measures the azimuthal angle from the wide side (indicated by W in the figure), to the narrow side (indicated by N in the figure). The variable $\hat {\xi }$ measures distance along the axis of the well. The thin dashed white line represents the interface between the two fluids, (with averaged concentration $\bar {c}=0.5$). Whereas the 2DGA model computes only the gap-averaged quantities only, the 3-D model computes variables across the annular gap. In order to compare with 2DGA results we post-process by averaging the solution variables from the 3-D model across the annular gap in the radial direction.

Figure 5. Gap-averaged concentration colour maps of case 4 (parameters: $e=0.6$; $Re=20$; $m=0.2$; $b=1000$) computed by (a) 2DGA model; (b) 3-D model.

From the comparison, both computations show a flat interface that evolves steadily along the annulus. Although there is a significant eccentricity the large buoyancy number leads to a steady displacement. In terms of predicting the overall behaviour of the (gap-averaged) 3-D code, the 2DGA is clearly effective. The main difference is in the intensity of the (red coloured) concentration close to the interface. The concentration is significantly less dispersed for the 2DGA model. Accordingly, the front of the displacing fluid of the 3-D model is slightly more advanced than that in figure 5(a). To the eye, it may appear that the right-hand sides of the snapshots of figure 5(b) have concentration that is significantly lower than figure 5(a), making the previous volumetric comparison in figure 4(a) appear doubtful. However, note that the differences in the colour maps are largest on the narrow side of the annulus, where the volume is smallest. Equally, there is a small amount of dispersion ahead of the $\bar {c}=0.5$ contour which is not easily seen. To quantify this dispersive aspect in an averaged way we plot a spatiotemporal map in figure 6. Concentrations from both models are averaged over the cross-section at fixed heights $\hat {\xi }$, to give $\bar {C}(\hat {\xi },\hat {t})$. The linear advance of the front is evident and very slightly faster in the 3-D case. Most of the dispersion is found just behind the front in the 3-D case.

Figure 6. Spatiotemporal colour maps of cross-section-averaged concentration $\bar {C}(\hat {\xi },\hat {t})$ of case 4 computed by (a) 2DGA model; (b) 3-D model.

If figure 5 were typical of all displacement flows, there would be little motivation for improving the 2DGA model. However, we shall see in other examples that although the 2DGA model captures the qualitative features of the 3-D flow, the quantitative differences can be significant. To understand the origin of these differences, we must first ask where the dispersion comes from, and secondly better understand the 2DGA model limitations. Considering the underlying assumptions of the 2DGA model, there are three principal simplifications: (i) neglect of inertial effects in the momentum balance; (ii) elimination of gradients of the viscous stress tensor, simplifying to a shear flow; (iii) simplification of the advective mass transport terms in (2.1). It is the last of these that has a clear dispersive effect.

Regarding the actual dispersion, note that in both models we have no diffusive terms in the concentration equations. For the ranges of typical parameters in these displacement flows, the Péclet numbers would typically be in the range $10^7\unicode{x2013}10^9$, which motivates the neglect of diffusion terms, although they can also be included (Maleki & Frigaard Reference Maleki and Frigaard2017). The numerical methods used here do in fact have well-controlled numerical diffusion, limiting significant smearing of the interface to a few mesh cells. However, intermediate concentrations that are computed at the interface are advected away strongly by secondary flows, regardless of diffusion. In order to effectively represent very small physical diffusivities would require a very fine mesh close to the interface, which is not practically feasible for general studies. Even with minimal diffusion, we would observe dispersion in gap-averaging the 3-D velocity profiles, which are Poiseuille-like.

Although the flow of figure 5(a) has a very sharp interface, illustrating the efficacy of the 2-D numerical scheme in limiting numerical diffusion, there is in fact considerable dispersion even here. To see this, note that with $e=0.6$, the far-field flow is significantly faster on the wide side of the annulus than the narrow side, but yet the displacement front advances steadily. As a rough guide, the ratio of wide to narrow side mean velocities in the far-field is $\approx (1+e)^2/(1-e)^2$ for Newtonian fluids. The steadily propagating front therefore requires a significant secondary azimuthal flow, moving fluid from wide to narrow (and vice versa). This secondary flow is localised close to the advancing front. Thus, considering figure 5(a,b), both models have secondary flows and consequent dispersion, but the effects are different.

Now we consider case 7, which has smaller eccentricity and buoyancy number. Figure 7 shows the gap-averaged concentration colour maps computed by both models. As the light blue colour ahead of the interface is not very clear to see, we plot yellow contour lines which represent $\bar {c}=0.01$. There are distinct differences between models. The 2DGA displacement is steady and piston-like at the front. For the 3-D model the interface ($\bar {c}=0.5$) is near-horizontal on the wide side, but begins to lag behind on the narrow side. Obvious dispersion is seen in the 3-D result. This is again primarily behind the displacement front, but we do see a region extending ahead of the front on the wide side. The dispersion is amplified by the small buoyancy number. Note that here $e=0.2$, but the effects of dispersion along the narrow side are far more significant than in case 4. There is no residual layer in the narrow side of the annulus, as the displacing fluids keep moving. There are also mild instabilities that start close to the interface in the narrow side, as the colour map shows at 85 s and 120 s, which do not apparently originate in the 2DGA physics, even though they are reminiscent of fingering-type instabilities.

Figure 7. Gap-averaged concentration colour maps of case 7 (parameters: $e=0.2$; $Re=20$; $m=0.5$; $b=10$) computed by (a) 2DGA model; (b) 3-D model.

To have a better understanding of the 3-D flow, figure 8 shows snapshots of the isosurface $c=0.2$, at different times for case 7, exploring the annulus section from 2.4 m to 3.6 m, i.e. concentrations of 20 % and above are enclosed within the isosurfaces shown. The colour map on the equally spaced cross-sections represents the dimensional axial velocity. We can clearly see the dispersion, along the centreline of each annular gap and faster on wide side. The dispersive velocity reached two times the imposed mean velocity. The gap-averaged velocity of the narrow side is smaller than that of the wide side, which results in the interface shown in figure 7(b), beginning to extend downstream. In addition, we show an annulus cross-section at $\hat {\xi }=3$ m (the middle of the annular sections of figure 8), to illustrate the azimuthal flow, see figure 9. The colour map shows the concentration distribution and black lines represent streamlines. We see a symmetric flow from narrow side to wide side before the main current passes (figure 9a,b), and the azimuthal flow is not that severe. Then, as the front displaces the secondary flow becomes stronger (figure 9c,d), and appears to lose its symmetry. The flow is mostly from wide side to narrow side, which appears to cause the instabilities observed in figure 7(b). As we gap-averaged the velocity field, the 2DGA model cannot represent the complex asymmetric patterns seen in figure 9(c,d), at least in its current form. The question is, whether such instabilities and asymmetries are caused by inertial or by gap-scale dispersive effects, which are not in the 2DGA model.

Figure 8. Snapshots of the 3-D displacement flow at different times for case 7, (parameters: $e=0.2$; $Re=20$; $m=0.5$; $b=10$): (a) $\hat {t}=55\,{\rm s}$; (b) $\hat {t}=65\,{\rm s}$; (c) $\hat {t}=75\,{\rm s}$; (d) $\hat {t}=85\,{\rm s}$. Colour map of equally spaced slices represents dimensional axial velocity. The wide side of the annulus is on the top.

Figure 9. Azimuthal flow of the interface located at $\hat {\xi }=3$ m at different times for case 7: (a) $\hat {t}=55\,{\rm s}$; (b) $\hat {t}=65\,{\rm s}$; (c) $\hat {t}=75\,{\rm s}$; (d) $\hat {t}=85\,{\rm s}$. The colour map shows the concentration distribution and yellow lines represent streamlines. The wide side of the annulus is on the top.

Further examples of the 2DGA model output and its comparison with 3-D simulations are given in the Appendix. As well as dispersion, these also show significant instabilities (more than those of figure 7), that are not captured by the 2DGA model. Together with the results above, these provide motivation for improving the 2DGA model, so as to better represent cementing displacements. The obvious missing ingredient is to account for dispersion due to the velocity profiles across the annular gap.

4. Accounting for gap-scale dispersion

To account for the gap-scale dispersion in the 2DGA we have two main options. First, we might resolve the flow on the scale of the gap using a simplified 3-D approach. For example, we might follow Tardy et al. (Reference Tardy, Flamant, Lac, Parry, Sri Sutama and Baggini Almagro2017) and Tardy (Reference Tardy2018) and discretise in the radial direction, while keeping the scaling arguments that lead to a shear flow at leading order. Evaluating the radial variation then allows us to advance the concentration using a velocity field that is 3-D and hence resolve the products of $c$ and the velocity directly. This leads to a model that is somewhere between 2-D and 3-D (say 2.5-D). The additional meshing in the radial direction means larger data sets to handle and slower computation, while still not resolving the inertial effects (as is done in a fully 3-D computation). Nevertheless, with care directed at the computational aspects, the 2.5-D approach of Tardy demonstrably works well for many cases.

The second option is to retain the 2-D description and to derive expressions that can express the gap average of the product of concentration and velocity, in a way that approximates the actual dispersive behaviour. The 2DGA model is compatible with the idea that dispersion is primarily in the direction of the gap-averaged flow. In Maleki & Frigaard (Reference Maleki and Frigaard2017) this idea is made more explicit in modelling the Taylor dispersion in turbulent annular flows. However, here the laminar flow is far from the Taylor dispersion regime. Instead, the flows are characterised as large Péclet number, small aspect ratio miscible displacement flows on the annulus gap-scale.

These regimes were studied by Yang & Yortsos (Reference Yang and Yortsos1997), who considered asymptotic solutions to the problem of miscible Newtonian displacement in long channels. Their leading-order transverse flow equilibrium (TFE) model combines a lubrication approximation for the velocity, with a convection–diffusion equation for the fluid concentration, containing the transverse diffusion only. Using this approach, a first-order conservation law was derived for the average concentration of displacing fluid ($\bar {c}$), along the plane channel,

(4.1)\begin{equation} \frac{\partial \bar{c}}{\partial t} + \frac{\partial }{\partial x}F_p = 0, \end{equation}

where $F_p$ is the flux of displacing fluid along the channel. In the limit of zero diffusion $F_p$ can be calculated explicitly and expressed as $F_p(\bar {c};m)$. Solutions of (4.1) exhibit a range of behaviours, including shocks and dispersive spreading. Similar flow regimes were analysed later by Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999) who extended the theory of Yang & Yortsos (Reference Yang and Yortsos1997) to include buoyancy effects. The propagation is governed again by (4.1), but now the flux is given by

(4.2)\begin{equation} F_p(\bar{c},m,U) = \frac{\bar{c}}{2}\frac{[(2m-3)\bar{c}^2+3]}{[1+(m-1)\bar{c}^3]}+\frac{\bar{c}^2(1-\bar{c})^3}{4U}\frac{[(4m-3)\bar{c}+3]}{[1+(m-1)\bar{c}^3]}. \end{equation}

The first term is identical to that of Yang & Yortsos (Reference Yang and Yortsos1997). Now $F_p$ depends also on additional parameter $U$, as follows:

(4.3)\begin{equation} U = \frac{3{\hat{\mu}_{1}\hat{w}_{0}}}{\hat{h}^{2}\hat{g}(\hat{\rho}_{1}-\hat{\rho}_{2})} .\end{equation}

In the model of Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999), $\hat {h}$ represents the channel half-width, and fluid 1 is displaced by fluid 2. The flow is vertical, in the direction of gravity and $U \ge 0$, i.e. density stable. In Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999) there is typographical error for (4.2): corrected in the PhD thesis (Lajeunesse Reference Lajeunesse1999), and above. The expression (4.2) has been analysed in Lajeunesse (Reference Lajeunesse1999), to predict flow regimes, which mostly agree with those observed experimentally. The question we now address is whether a similar methodology can be effectively used in our 2-D flows, where the channel width and buoyancy direction varies, and where the fluids may also be non-Newtonian.

4.1. The Hele-Shaw setting

If one considers primary cementing in a vertical casing to represent a Hele-Shaw cell of varying gap, it is tempting to resolve the displacement flow in vertical slices along the length of the well. On each such slice we might expect the above analysis to be valid and deliver a good approximation of the gap-scale dynamics. While correct, this neglects the effects of azimuthal secondary flows which can be significant, as we have seen; fortunately, there is a better solution.

In the derivation of the 2DGA model, (2.1)–(2.3), the Navier–Stokes equations are made dimensionless using a different length scale in the radial direction compared with that in azimuthal and axial directions. This leads to an asymptotic approximation (in terms of the aspect ratio), which to leading order is a shear flow, i.e.

(4.4)$$\begin{gather} 0 ={-} \frac{\partial p}{\partial y}, \end{gather}$$

(4.5)$$\begin{gather}0 ={-} \frac{1}{r_a}\frac{\partial p}{\partial \phi} +\frac{\partial}{\partial y} \tau_{\phi y} + \frac{\rho \sin \beta \sin {\rm \pi}\phi}{{F^2}} , \end{gather}$$

(4.6)$$\begin{gather}0 ={-} \frac{\partial p}{\partial \xi } + \frac{\partial}{\partial y} \tau_{\xi y} - \frac{ \rho \cos \beta }{{F^2}}, \end{gather}$$

where $y \in [-H(\phi,\xi ),H(\phi,\xi )]$ is the scaled radial (gap-)coordinate. These equations have scaled the pressure gradients in order to balance with the leading-order shear flow. The density above has been scaled with the density of fluid 1. We deviate slightly from Maleki & Frigaard (Reference Maleki and Frigaard2017) above, since we have chosen to not subtract the static pressure field of fluid 1 from the pressure before scaling. This is because evaluating buoyancy effects requires more care. The leading-order strain rate is

(4.7)\begin{equation} \dot{\gamma} \sim \left[ \left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial w}{\partial y}\right)^2\right]^{1/2}. \end{equation}

Thus, to leading order the constitutive laws depend only on the $y$-derivatives of both $v$ and $w$, and on the fluid present. By integrating across the gap, it is found that the motion is driven in the direction of the modified pressure gradient,

(4.8)\begin{equation} \boldsymbol{G} = (G_{\phi},G_{\xi}) \equiv \left(- \frac{1}{r_a}\frac{\partial p}{\partial \phi} + \frac{ \rho \sin \beta \sin {\rm \pi}\phi}{{F^2}}, -\frac{\partial p}{\partial \xi } - \frac{ \rho \cos \beta }{{F^2}} \right). \end{equation}

In other words, the vector of averaged velocities $(\bar {v},\bar {w})$ is parallel in the $(\phi,\xi )$ plane to $\boldsymbol {G}$.

The difficulty comes in resolving the averaged velocities. In Bittleston et al. (Reference Bittleston, Ferguson and Frigaard2002) and Maleki & Frigaard (Reference Maleki and Frigaard2017), the simplification was made that the concentrations were constant across the annular gap. With this assumption, the density $\rho$ and any rheological properties should be defined as functions of the (mean) fluid concentrations. For the density a simple volume average is appropriate. Having specified the rheology of the mixture, the relevant closure expression can also be derived, based on a plane Poiseuille flow. However, it is clear that these simplifications negate dispersion that will occur due to layering of the fluids across the gap during displacement. Equally, buoyancy effects are limited to gradients of the mean density and not individual fluids.

4.2. Dispersive Hele-Shaw setting

In principle, the same method as Bittleston et al. (Reference Bittleston, Ferguson and Frigaard2002) and Maleki & Frigaard (Reference Maleki and Frigaard2017) can be employed here, but with different assumptions made regarding the distribution of fluids across the gap. For simplicity we assume only two fluids: fluid 1 displaced by fluid 2. We assume that the displacing fluid can disperse (symmetrically) down the centre of the local annular gap. Given the large $Pe$ we assume that fluid 2 occupies the centre of the channel, between interfaces at $y = \mp y_i$; fluid 1 occupies the wall layers: $[-H,-y_i)$ and $(y_i,H]$ (see figure 10). The density $\rho$ in (4.5) and (4.6) will now be constant in each fluid layer, again allowing integration, but now with a jump in fluid properties at $y = \mp y_i$.

Figure 10. Schematic of dispersive Hele-Shaw model setting: (a) narrow annulus, with slice along the local flow direction $\boldsymbol {s}$; (b) assumed configuration of fluid 1 (blue) and fluid 2 (red); local concentration $c = y_i/H$.

The concentration of displacing fluid 2 in the gap we shall denote $\bar {c}(\phi,\xi,t)$, which is given geometrically by $\bar {c} = y_i/H$. The evolution of $\bar {c}$ (neglecting diffusive terms), is governed by

(4.9)\begin{equation} \frac{\partial}{\partial t}[Hr_{a}\bar{c}]+ \boldsymbol{\nabla}_a \boldsymbol{\cdot} \boldsymbol{q} = 0, \end{equation}

where the flux functions are

(4.10)\begin{equation} \boldsymbol{q} = (q_\phi , q_\xi ) = r_a\left( \int^{y_i}_0 v\,{{\rm d} y}, \int^{y_i}_0 w\,{{\rm d} y} \right). \end{equation}

Note that if $(v,w)$ are approximated by their gap-averaged quantities, then (4.9) is exactly (2.1). The questions therefore are, how to evaluate the above fluxes and are we able to specify them in terms of gap-averaged quantities?

The shear stresses in (4.5) and (4.6) are coupled only through the effective viscosity of each fluid, $\eta _1$ and $\eta _2$, which depend on $\dot {\gamma }$. On assuming symmetry about $y=0$, integrating (4.5) and (4.6) for $y \in [0,H]$, and imposing continuity of the shear stresses at $y=y_i$, we find

(4.11)\begin{gather} \frac{\partial v}{\partial y} = \left\{ \begin{array}{@{}ll} \displaystyle{\dfrac{y}{\eta_2(y)}\left( \dfrac{1}{r_a}\dfrac{\partial p}{\partial \phi} - \dfrac{\rho \,f_{\xi}}{r_a} \right) + \dfrac{y}{\eta_2(y)} \dfrac{(H-y_i) \Delta \rho \,f_{\xi} }{H r_a} } , & y \in [0,y_i],\\ \displaystyle{ \dfrac{y}{\eta_1(y)} \left( \dfrac{1}{r_a}\dfrac{\partial p}{\partial \phi} - \dfrac{\rho \,f_{\xi}}{r_a} \right) + \dfrac{(H-y)}{\eta_1(y)} \dfrac{y_i\Delta \rho \,f_{\xi} }{Hr_a} }, & y \in (y_i,H], \end{array} \right. \end{gather}

(4.12)\begin{gather} \frac{\partial w}{\partial y} = \left\{ \begin{array}{@{}ll} \displaystyle{\dfrac{y}{\eta_2(y)}\left(\dfrac{\partial p}{\partial \xi} + \dfrac{\rho \,f_{\phi}}{r_a} \right) - \dfrac{y}{\eta_2(y)} \dfrac{(H-y_i) \Delta \rho \,f_{\phi} }{H r_a} } , & y \in [0,y_i],\\ \displaystyle{ \dfrac{y}{\eta_1(y)} \left(\dfrac{\partial p}{\partial \xi} + \dfrac{\rho \,f_{\phi}}{r_a} \right) - \dfrac{(H-y)}{\eta_1(y)} \dfrac{y_i\Delta \rho \,f_{\phi} }{H r_a} }, & y \in (y_i,H], \end{array} \right. \end{gather}

where $\Delta \rho = \rho _1-\rho _2$ is the density difference, $\rho = (1-y_i/H) \rho _1+ (y_i/H)\rho _2$ is the mean density and

(4.13)\begin{equation} \boldsymbol{f} = (\,f_{\phi},\,f_{\xi}) = \left( \frac{r_a \cos \beta }{{F^2}},\frac{r_a \sin \beta \sin {\rm \pi}\phi}{{F^2}}\right) . \end{equation}

Observe that the velocity gradients depend on two terms: a pressure gradient that is modified with the mean static pressure gradient, cf. (4.8), and then on a buoyancy term that scales with $\Delta \rho$. Note also that $\rho _1 = 1$, due to scaling with the density of fluid 1.

These expressions for the velocity gradients can be integrated further, outwards from the wall at $y=H$, where the velocity is zero, to give $v(y)$ and $w(y)$. A further integration over $[0,H]$ produces expressions for $(H\bar {v},H\bar {w})$. The point of this laborious exercise is to see that the integral expressions that multiply the modified pressure gradient terms and those that multiply the buoyancy terms, are identical in the two directions,

(4.14)\begin{equation} (H\bar{v},H\bar{w}) ={-}\mathcal{I}_1 \left(\frac{1}{r_a}\frac{\partial p}{\partial \phi}- \frac{\rho \,f_{\xi}}{r_a} , \frac{\partial p}{\partial \xi} + \frac{\rho \,f_{\phi} }{r_a} \right) + \mathcal{I}_2 \frac{\Delta \rho }{Hr_a}(\,{-}f_{\xi}, \,f_{\phi} ). \end{equation}

This is valid for any generalised Newtonian fluid. The expressions for $\mathcal {I}_1$ and $\mathcal {I}_2$ are defined as

(4.15)$$\begin{gather} \mathcal{I}_1= \left[ y_i \int^H_{y_i} \frac{\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y} + \int_0^{y_i} \int^{y_i}_y \frac{\tilde{y}}{\eta_2(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d}y} + \int_{y_i}^H \int^H_y \frac{\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d} y} \right], \end{gather}$$

(4.16)$$\begin{gather}\mathcal{I}_2 = \left[ (H-y_i) \int_0^{y_i} \int^{y_i}_y \frac{\tilde{y}}{\eta_2(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d} y} + y_i^2 \int^H_{y_i} \frac{H-\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y} + y_i \int_{y_i}^H \int^H_y \frac{H-\tilde{y}}{\eta_1(\tilde{y})}\,{\rm d}\tilde{y}\,{{\rm d} y} \right]. \end{gather}$$

These expressions depend on the constitutive laws for both fluids, as well as $y_i$ and $H$. For most generalised Newtonian fluids we will not be able to evaluate $\mathcal {I}_1$ and $\mathcal {I}_2$, except numerically.

Considering (4.14), in the absence of the buoyancy term, e.g. if $\Delta \rho = 0$, we can see that the mean flow follows the direction of the modified pressure gradient, as in Bittleston et al. (Reference Bittleston, Ferguson and Frigaard2002), and see (4.8). The mean density is again volumetrically averaged. With a density difference, however, the flow direction is influenced independently by the buoyancy vector: the second term on the right-hand side of (4.14). The form of (4.14) has the implication that the shear flow is a 1-D flow, in the direction of the gap-averaged velocity. Although we do not know this direction a priori, we may proceed assuming that this direction is known, solve the 1-D flow, then effectively calculate $\mathcal {I}_1$ and $\mathcal {I}_2$ to determine the flow direction from (4.14).

4.2.1. The 1-D shear flow

We assume that $(\bar {v},\bar {w})$ flows locally in direction $\boldsymbol {s} = (s_\phi,s_\xi )$. Taking the dot product of (4.5) and (4.6) with $\boldsymbol {s}$ we have the following 1-D momentum balance:

(4.17)\begin{equation} \frac{\partial }{\partial y} \tau_{sy} = \left\{ \begin{array}{@{}ll} \displaystyle{ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) - \dfrac{(H-y_i) \Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) }, & y \in [0,y_i],\\ \displaystyle{ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) + \dfrac{y_i \Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) }, & y \in (y_i,H]. \end{array} \right. \end{equation}

This type of buoyant two-layer flow model can generally be resolved numerically in an efficient way for two generalised Newtonian fluids. For example, see Zare, Roustaei & Frigaard (Reference Zare, Roustaei and Frigaard2017) for Bingham fluids, or Taghavi et al. (Reference Taghavi, Seon, Martinez and Frigaard2009) for a different two-layer model for Herschel–Bulkley fluids. The problems are tractable due to various monotonicity results. Here, however, we focus on the case of Newtonian fluids for simplicity.

We first integrate to find the shear stress, which varies linearly in each layer,

(4.18)\begin{equation} \tau_{sy} = \left\{ \begin{array}{@{}ll} \displaystyle{\left[ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right]y - (H-y_i)y \dfrac{\Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-f_{\xi}s_{\phi}) }, & y \in [0,y_i],\\ \displaystyle{ \left[ \dfrac{\partial p}{\partial s} + \dfrac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right] y - (H-y)y_i \dfrac{\Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) }, & y \in (y_i,H] . \end{array} \right. \end{equation}

Note that we are not interested (yet) in the actual velocity solution $u_s(y)$, but in the areal flux $H\bar {u}_s$. This can be evaluated for Newtonian fluids by integrating by parts,

(4.19)\begin{align} H\bar{u}_s &={-} \int_0^{y_i} y \frac{\partial u_s}{\partial y}\,{{\rm d} y} - \int_{y_i}^H y \frac{\partial u_s}{\partial y}\,{{\rm d} y} ={-} \int_0^{y_i} y \frac{\tau_{sy}}{\eta_2}\,{{\rm d} y} - \int_{y_i}^H y \frac{\tau_{sy}}{\eta_1}\,{{\rm d} y} \nonumber\\ & ={-} \left( \frac{y_i^3}{3 \eta_2 } + \frac{H^3 - y_i^3}{3 \eta_1 } \right) \left[ \frac{\partial p}{\partial s} + \frac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right] \nonumber\\ &\quad + \frac{\Delta \rho }{H r_a} (\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi})\left[ \frac{(H-y_i)y_i^3}{3 \eta_2 } + \frac{y_i(H-y_i)^2(H+2 y_i)}{6 \eta_1} \right]. \end{align}

We may also evaluate $H\bar {u}_s$ from (4.14),

(4.20)\begin{align} H\bar{u}_s &= \boldsymbol{s} \boldsymbol{\cdot} (H\bar{v},H\bar{w}) \nonumber\\ & ={-} \mathcal{I}_1 \left[ \frac{\partial p}{\partial s} + \frac{\rho}{r_a}(\,f_{\phi}s_{\xi}-\,f_{\xi}s_{\phi}) \right] + \mathcal{I}_2 \frac{\Delta \rho }{Hr_a}(\,f_{\phi} s_{\xi} -\,f_{\xi} s_{\phi}) . \end{align}

Comparing these two expressions, we see that for a Newtonian fluid pair,

(4.21a,b)\begin{equation} \mathcal{I}_1 = \left(\frac{y_i^3}{3 \eta_2 } + \frac{H^3 - y_i^3}{3 \eta_1 } \right),\quad \mathcal{I}_2 = \left[ \frac{(H-y_i)y_i^3}{3 \eta_2 } + \frac{y_i(H-y_i)^2(H+2 y_i)}{6 \eta_1} \right] . \end{equation}

Power law and yield-stress fluid pairs will be dealt with in a subsequent study, requiring computation. Having solved for the 1-D shear flow, we may now go back to (4.14) and derive the modified 2-D model.

4.2.2. The D2DGA model

The gap averaged stream function is defined as ${\nabla }_a \varPsi = (2H\bar {w},-2H \bar {v})$, which can be substituted into (4.14). To derive the D2DGA model, we rearrange (4.14) to give two equations for the two pressure gradient components. We then cross-differentiate to eliminate the pressure,

(4.22)\begin{align} 0 &= \boldsymbol{\nabla}_a \boldsymbol{\cdot} \left[ \left( \frac{r_a}{2\mathcal{I}_1} \right){\nabla}_a \varPsi \right] + \boldsymbol{\nabla}_a \boldsymbol{\cdot}\left[ \left( \rho - \frac{\Delta \rho}{H}\frac{\mathcal{I}_2}{\mathcal{I}_1} \right) \boldsymbol{f} \right] \nonumber\\ & = \boldsymbol{\nabla}_a \boldsymbol{\cdot} [\boldsymbol{S} + \boldsymbol{b}] . \end{align}

We observe that (4.22) is an elliptic Poisson problem, of the same form as in Bittleston et al. (Reference Bittleston, Ferguson and Frigaard2002) and Maleki & Frigaard (Reference Maleki and Frigaard2017). However, the assumption of a layered flow gives rise to significant differences.

Considering first $\boldsymbol {S}$, we see that for a Newtonian fluid,

(4.23)\begin{equation} \boldsymbol{S} = \left( \frac{r_a}{2\mathcal{I}_1} \right){\nabla}_a \varPsi = \frac{3r_a }{2H^3} \frac{{\nabla}_a \varPsi}{\displaystyle{\frac{\bar{c}^3}{\eta_2} + \frac{(1-\bar{c}^3)}{\eta_1}}},\quad \Rightarrow\quad \tau_w = \frac{3}{2H^2} \frac{|{\nabla}_a \varPsi |}{\displaystyle{\frac{\bar{c}^3}{\eta_2} + \frac{(1-\bar{c}^3)}{\eta_1}}} , \end{equation}

as opposed to $\boldsymbol {S} = [3 r_a\eta /2H^3] {\nabla }_a \varPsi$ for the 2DGA models. If the 2DGA approach were to yield the same $\boldsymbol {S}$, we see that a specific form of viscosity interpolation would be needed. In practice a linear interpolation has been used for 2DGA computations, with the assumption that the fluids are mixed.

Turning now to the buoyancy vector $\boldsymbol {b}$, we see that there are two parts in (4.22). Firstly, if we write $\rho = (1-\bar {c}) \rho _1 + \bar {c} \rho _2$, this helps us to see that the first part of $\boldsymbol {\nabla }_a \boldsymbol {\cdot} \boldsymbol {b}$ represents changes in the mean density in the direction $\boldsymbol {f}$. These buoyancy gradients are exactly those considered in the models of Bittleston et al. (Reference Bittleston, Ferguson and Frigaard2002) and Maleki & Frigaard (Reference Maleki and Frigaard2017). The second part results specifically from the layered flow. If we write $\rho = \rho _2 + (1-\bar {c}) \Delta \rho$, we see that the buoyancy gradients from the mean density term have similar order of magnitude as those that result from the second part. In other words, the effect of the layered flow on buoyancy is significant. In more detail, we have

(4.24)\begin{equation} \boldsymbol{\nabla}_a \boldsymbol{\cdot} \boldsymbol{b} = \boldsymbol{\nabla}_a \boldsymbol{\cdot} \left[ \rho_2 + \Delta \rho \left( (1-\bar{c}) - \left[ \frac{m\bar{c}^3(1-\bar{c}) + \bar{c}(1-\bar{c})^2(0.5+\bar{c})}{m\bar{c}^3 + 1 - \bar{c}^3} \right] \right) \right] \boldsymbol{f} . \end{equation}

The gradients of $\boldsymbol {f}$ will typically be small as changes in the well orientation occur slowly.

With regard to transport of the fluid concentrations, we need to evaluate the fluxes. For Newtonian fluids this is straightforward, although algebraically long. The pressure gradients are eliminated using (4.14) and there is some manipulation using the expressions $\mathcal {I}_1$ and $\mathcal {I}_2$. This simplifies to

(4.25)$$\begin{gather} (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} \left[ \frac{m\bar{c}^2 + 1.5(1-\bar{c}^2)}{m\bar{c}^3 + 1-\bar{c}^3}\right] + \frac{\Delta \rho H^3}{6 \eta_2} \mathcal{I}_3(\bar{c},m)[{-}f_\xi,f_\phi ], \end{gather}$$

(4.26)$$\begin{gather}\mathcal{I}_3(\bar{c},m) = \frac{\bar{c}^2(1-\bar{c})^3[4m\bar{c} + 3(1-\bar{c})]}{2m[m\bar{c}^3 + 1 - \bar{c}^3 ]}. \end{gather}$$

We observe that the first term is very similar to the model of Maleki & Frigaard (Reference Maleki and Frigaard2017), except that the transport term in this model is

(4.27)\begin{equation} \frac{r_a}{2} \left(- \frac{\partial \varPsi}{\partial \xi } , \frac{1}{r_a}\frac{\partial \varPsi}{\partial \phi } \right) \bar{c} . \end{equation}

We see that this is amplified by the concentration dependent flux term,

(4.28)\begin{equation} \frac{m\bar{c}^2 + 1.5(1-\bar{c}^2)}{m\bar{c}^3 + 1-\bar{c}^3}. \end{equation}

On inspection, we see that this amplification factor is the first term in (4.2), as derived by Yang & Yortsos (Reference Yang and Yortsos1997). Unlike these earlier studies, in the annulus the gap width and local mean velocity may now vary.

The second term in (4.25) is also related to (4.2). If for example we simplify to have a single component of $\boldsymbol {f}$, e.g. in a vertical well with $\beta = 0$, then we see that $f_\xi =0$. The $\phi$-component of $\boldsymbol {q}$ consists of only the first dispersive term of (4.25), whereas the $\xi$-component is augmented by buoyancy. The $\bar {c}$-dependence of $\mathcal {I}_3(\bar {c},m)$ mimics that of the second term in (4.2) and the expressions are identical if

(4.29)\begin{equation} U = \frac{3m \eta_2 {F^2}}{\Delta \rho H^3 r_a}. \end{equation}

If the annulus is uniform and concentric ($H=r_a=1$), we find that these terms are identical.

6. Summary

This paper has studied primary cementing displacement flows in vertical eccentric annuli. These flows are challenging to compute, due to the excessively long length scales of typical wells (${\sim} 10^3$ m), compared with the annular gap (${\sim} 10^{-2}$ m). Ensuring a mesh size needed for accuracy, both in representing the interface and on the scale of the annular gap, while also integrating in time long enough for the displacement front to travel the length of the well, makes 3-D non-Newtonian computations unfeasible at present. Nevertheless, it is important to understand the effectiveness of the full displacement process for engineering design purposes. By design we mean selection of the annulus eccentricity, the fluid densities and rheologies, the volumes and flow rates pumped. Each of these may only be varied within operational constraints. Because of this need, 2DGA models have become popular over the past 20 years.

In this paper we have used 3-D transient computations in annuli of length equal to ${\sim}10^3$ annular gap widths, to represent the full complexity of the displacement flow, and compare these with 2DGA results. While 2DGA models can predict many qualitative phenomena well, (e.g. steady/unsteady, stable/unstable, residual narrow side channel or not), there are also quantitative discrepancies in most cases. Computed 3-D simulations predict fluid concentrations that are clearly dispersive and that result in lower displacement efficiencies than 2DGA flows.

To account for gap-scale dispersion in the 2DGA framework, we rederive the model, starting from the reduced shear flow equations. We assume that the displacing fluid can symmetrically displace down the centre of the local annular gap, instead of assuming a fully mixed concentration across the gap (Bittleston et al. Reference Bittleston, Ferguson and Frigaard2002; Maleki & Frigaard Reference Maleki and Frigaard2017). This leads to a buoyant two-layer flow model for the 1-D shear flow. The D2DGA model is obtained after solving the shear flow equations and eliminating the pressure components. This model derivation is presented for any two generalised Newtonian fluids. However, for many of these a numerical solution on the gap-scale would be needed. On the other hand, for the two Newtonian fluids that are used in our examples, the model can be found analytically.

The D2DGA equation for the stream function is an elliptic Poisson problem, just as in the 2DGA model. However, there are differences in closure expressions for both the wall shear stress and buoyancy. With regard to the mass transport equation, we obtain a newly defined gap-averaged flux that can approximate the actual dispersive behaviour. An additional term representing buoyancy effects in the layered flow also arises, compared with the 2DGA model. As the underlying mathematical structure of the 2DGA and D2DGA models is the same, computational times are similar, with some additional cost from evaluating the new closure functions.

Of course, the D2DGA model is not the first study of dispersive effects in duct-like displacements. In the classical study of Saffman & Taylor (Reference Saffman and Taylor1958) the authors explore viscosity (mobility) ratio effects on instabilities in immiscible Hele-Shaw displacements: viscous fingering. Much of this paper concerns viscous fingering in the case in which there is a single fluid present in the through-thickness direction. However, the Appendix contains a modified analysis of the situation with incomplete displacement, i.e. the displacing fluid has a multilayer structure. As here, this renders the mean velocity concentration dependent. Secondly, multiphase porous media flow models applied in oil recovery often employ flux functions that describe the wetting of different phases, dating from Buckley & Leverett (Reference Buckley and Leverett1942). These expressions are either empirically derived or rely on simplified models, such as the Hele-Shaw cell or tube bundles, in order to develop mathematical closure expressions. The wetting closures are then applied to rather complex actual pore geometries. Thus, our approach is certainly related to these larger application areas.

Perhaps we are fortunate in that the narrow eccentric annulus of primary cementing forms a prototype that has systematic and understandable geometry variation (effectively the permeability), allowing for the modelling approach to be effective. As displacement flows in cementing are also (predominantly) miscible and large Péclet number, our extension of Yang & Yortsos (Reference Yang and Yortsos1997) and Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999) seems the appropriate modelling approach. For uniform, vertical concentric annuli, the new closure expressions we derive are identical with those of Yang & Yortsos (Reference Yang and Yortsos1997) and Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999), who considered high Péclet number miscible Newtonian displacement flows in long vertical channels.

Through detailed comparisons between the D2DGA model and the 3-D model, we find significant improvement of the D2DGA model in capturing gap-scale dispersion. Moreover, some instabilities are predicted by D2DGA model, similar to those in the 3-D model. The D2DGA model remains computationally efficient and appears able to represent the different scales of dispersion, consistent with (averaging of) 3-D simulation results. The one deficiency is in modelling density-unstable displacement flows. Although the D2DGA prediction may be an improvement, quantitative prediction remains poor, as might be expected.

Finally, the effects of the dimensionless parameters have been explored, in terms of breakthrough times and displacement efficiency. The buoyancy number is the key factor in vertical primary cementing, once the eccentricity is fixed. Negative buoyancy number should be strictly avoided to prevent highly ineffective cementing and compromised well integrity. A large enough buoyancy number can be effective for Newtonian fluids, even if the annulus is highly eccentric. It helps with displacing the fluids on the narrow side and reducing residual layers caused by eccentricity. Small viscosity ratio (displaced/displacing) is also advisable. In the D2DGA model this directly decreases the dispersion between displacing and displaced fluids, essentially by acting on the drainage layers adjacent to the walls. Especially when the buoyancy number is small, the viscosity ratio affects dispersion more significantly, i.e. large viscosity ratio $m$ will lead to higher dispersion and less efficiency. Eccentricity affects the secondary flow and thus contributes to uniform/non-uniform distribution in the azimuthal direction. Large eccentricity is not beneficial for complete and efficient displacement. Here we only studied Newtonian fluids, but for yield stress fluids it is also well known that these fluids can become stuck on the narrow side for large $e$, (McLean et al. Reference McLean, Manry and Whitaker1967; Couturier et al. Reference Couturier, Guillot, Hendriks and Callet1990). As for the effect of Reynolds number, we did not find that it plays an essential role and especially with large $b$. Although operationally speaking, this simply affirms an old story (Couturier et al. Reference Couturier, Guillot, Hendriks and Callet1990; Frigaard & Pelipenko Reference Frigaard and Pelipenko2003; Pelipenko & Frigaard Reference Pelipenko and Frigaard2004c; Nelson & Guillot Reference Nelson and Guillot2006), this paper gives far greater detail. The real novelty is in being able to access predictions of the dispersive behaviour that agree well with the more expensive 3-D models.

The next step in our study is to consider how effective the D2DGA model is in predicting experimental results from our flow loop. Thereafter, the model will be applied to inclined and horizontal annuli, which can be achieved easily by changing the direction of the buoyancy force ($\beta \neq 0$). Here we already have the recent studies of Renteria & Frigaard (Reference Renteria and Frigaard2020) and Sarmadi et al. (Reference Sarmadi, Renteria and Frigaard2021) to compare with. The next step is to develop the model targeting non-Newtonian fluids, as are more common in field situations with laminar flows. Although a more complicated computation is needed to resolve the gap-scale flow, e.g. for Herschel–Bulkley fluids, the basic principle is applicable as outlined in this paper. Another interesting question is whether the new buoyancy, wall shear stress and flux closures can be analysed in a similar way to Pelipenko & Frigaard (Reference Pelipenko and Frigaard2004c) and Maleki & Frigaard (Reference Maleki and Frigaard2020), to rapidly classify cementing displacements in a way that is suitable for machine learning approaches.