3.1.1. Classifying steady and unsteady advective behaviour via a similarity transform

The first task is meaningful to classify displacements as either a steady front or unsteady. Although ideally this should mean that the front advances at the mean pumping speed all around the annulus, dispersion will lead to a faster arrival of the displacing fluid. Equally, we have three methods for assessing the flow: (i) an imperfect experiment; (ii) 3-D simulation; and (iii) D2DGA simulations. Our approach is as follows. For the experiments, we begin by selecting images taken from the frontal view of the first and last sections at different time steps. We average the image $\bar {C}_x(\hat {y},\hat {z},\hat {t})$ with respect to $\hat {y}$ to give $\bar {C}(\hat {z},\hat {t})$. Since we already averaged across the annular cross-section in the $x$-direction, the quantity $\bar {C}(\hat {z},\hat {t})$ represents the (experimental) mean concentration of displacing fluid at each distance, obtained by a cross-sectional area average. For the 3-D results, we have already averaged first in $\hat {x}$ to give $\bar {c}_x$, (comparable to the experimental image), and now also average with respect to $\hat {y}$ to give the cross-sectional area average $\bar {c}(\hat {z},\hat {t})$. Lastly, we average the D2DGA results with respect to the azimuthal direction. We also denote the latter as $\bar {c}(\hat {z},\hat {t})$. Since the radial coordinate was already averaged to predict $\bar {c}_r$, it follows that $\bar {c}(\hat {z},\hat {t})$ is also a cross-sectional area average.

Terms such as ‘front’ and ‘interface’ are always difficult when dealing with miscible displacements, as formally there is no interface. Equally, when we have significant dispersion, even for high-Péclet-number flows, the location of a nominal interface can become progressively smeared. Our visualization presents an image that is itself a type of average, hence vulnerable to dispersive smearing. A legitimate question is whether we can define a displacement front meaningfully at all that is applicable to the range of observed behaviours. To address this, we have focused on extracting the advective character of the flows. By plotting $\bar {C}(\hat {z},\hat {t})$ or $\bar {c}(\hat {z},\hat {t})$ against the variable $\hat {z}/\hat {t}$ for successive images (time steps), we are able to discern if the displacement is primarily advective. In the case where the flow is advection dominated, the mean concentrations will collapse onto a master curve in $\hat {z}/\hat {t}$. Denoting the normalized similarity variable by $w_f=(\hat {z}/\hat {t})/\hat {w}_{0}$, we find the master curve for $\bar {C}(w_f)$ or $\bar {c}(w_f)$ (3-D and D2DGA). The values of $w_f$, plotted as a function of cross-sectionally averaged concentration, show the range of front velocities found via each method.

This method is equally applicable to a nice clean interface or a very diffuse one, as we will see below. When the similarity variable is defined for all $\bar {C}$ (or $\bar {c}$), then the displacement front can be considered as the union of all the wavespeeds $w_f$. If the actual images were visualized, this visual representation might appear to be steady or otherwise, relatively clean or smeared. The ability to handle these different behaviours is the main attraction of using this method. The disadvantage is some ambiguity in defining the displacement front via $w_f$, as for cases of large dispersion, the actual image may not appear front-like at all.

Typically, the collapse of the averaged concentration onto the similarity variable $\hat {z}/\hat {t}$ is imperfect at both small and large mean concentrations, requiring filtering. As a threshold for the experimental data, we adopt the same criterion as Renteria & Frigaard (Reference Renteria and Frigaard2020), that is, we categorize flows as having a steady front by

(3.1)\begin{equation} {\rm \Delta} w_f = w_f(\bar{C} =0.3)-w_f(\bar{C} =0.7) \le 0.1. \end{equation}

The flows are classified as unsteady otherwise. When steady, we refer to the steadily moving segment between $\bar {C}=0.3$ and $\bar {C}=0.7$ as the main front. The front may of course extend beyond these threshold values, which are conservatively chosen. We start with a case of zero buoyancy force (EXP 5) for which $(e, b, m, Re) = (0, 0, 0.2, 50)$. Figure 4 shows a comparison of experiments, depth-averaged 3-D and D2DGA results. Figure 4(a–c) present the plotting of $\bar {C}$, $\bar {c}$ and $\bar {c}$ versus the similarity variable $w_f$, at a series of time steps, for the three methods. The experimental profiles are obtained only from the first and last sections, while the computational ones are computed from the entire annulus. For clarity, the mean concentrations when the front enters the last section of the annulus are plotted in green for the experiment, typically representing the most converged profiles, and the last profile computed is plotted for the simulations in green. Figure 4(d–f) display example displacement images from the first and last sections at two distinct time steps. The steps are chosen to be illustrative, and represent the same times for experiments and models. The same figure format is adopted subsequently.

Figure 4.

Case EXP 5 with $(e, b, m, Re) = (0, 0, 0.2, 50)$. Panels (a–c) plot $\bar {C}$, $\bar {c}$ and $\bar {c}$, from experiment, 3-D simulation and D2DGA models, respectively. The curves are collapsed onto the similarity variable $w_f=(\hat {z}/\hat {t})/\hat {w}_{0}$, representing the scaled front speed. The black and green curves in (a) represent the profiles of averaged concentration from the first and last section of the annulus from experimental results. The black profiles in (b,c) are computed from the entire annulus. Panels (d–f) show illustrative images ($\bar {C}_{x}$) from the frontal view of the experiment in (d), the depth averaged 3-D computation $\bar {c}_{x}$, in (e) and the D2DGA radially averaged concentrations $\bar {c}_r$ in (f). The times for experiment and simulations are the same. One image is taken from when the front is in the first section of the annulus and the other one from the fourth section of the annulus.

We observe that for the 3-D and D2DGA computational results, the data collapse very well onto a master curve $\bar {c}(w_f)$. However, the experimental data do not collapse well for this particular experiment (figure 4a), instead showing pronounced fluctuations. Additionally, the experimental curves appear to show larger maximum front velocity, whereas the models both have a peak front speed of approximately 1.5. The latter is consistent with dispersion along a plane channel. It is apparent that for this iso-dense case, ${\rm \Delta} w_f = w_f(\bar {C}=0.3)-w_f(\bar {C}=0.7)> 0.1$ from experimental and computational results, which indicates a typical unsteady flow.

Why the experimental results do not collapse onto a master curve $\bar {C}(w_f)$ and why we see peak front velocities significantly larger is probably due to a combination of factors. Reflections have a significant impact on the pixel values recorded. Therefore, the experimental results are not as distinct as the computational results in any of our experiments. There are also likely to be minor non-uniformities in the apparatus, deviating from truly concentric. Although the displacing fluid is more viscous ($m = 0.2$), this is the only stabilizing feature of the flow that might lead to a uniform and efficient displacement.

Looking at figure 4(d,e), we observe darker profiles close to the edges of the image, which are effectively the displacement front viewed in the annular gap on each side. This is consistent with a Poiseuille-type flow on the gap scale: with a dispersive spike that advances along the centre of the annular gap, leaving behind residual layers on both walls. The depth-averaged 3-D results capture this behaviour as well as the experiments. In the D2DGA model, this effect is implicitly included in the model. Whereas the thin spike produces near identical behaviour in the two models, in the experiment, we see that the spike is intermittent on the two sides and the images are also not uniform in $\hat {y}$, with some channelling occurring in the last section; see the clear white zone at the top left of the last section of figure 4(d). This asymmetry likely allows the (averaged) front to advance faster than in either model.

Now if we turn to look at a case (EXP 96) with large buoyancy number ($b=750$), figure 5 shows a quite different result. First, we can observe a flat front and complete displacement from the first to last sections through all of the results shown in figure 5(d–f). Second, the average concentration in the $\hat {y}$-direction collapses onto a master curve $\bar {c}(w_f)$, which clearly shows a steadily advancing front (with speed close to $w_f = 1$); see figure 5(a–c). This indicates an ideal steady front displacement with minimal dispersion. For the experimental results in figure 5(a), the series of curves from the first section has not converged as completely as those from the last section. For the model data, the last curve $\bar {c}(w_f)$ plotted has been coloured green, from which we can see the convergence. In all cases, the criterion (3.1) is satisfied. Compared with the previous example, it is evident that the buoyancy force is primarily responsible for limiting gap-scale dispersion, leading to the steady front displacement.

Figure 5.

Case EXP 96 with $(e, b, m, Re) = (0, 750, 0.2, {29})$. Panel (a–c) plot $\bar {C}$, $\bar {c}$ and $\bar {c}$ from experiment, 3-D simulation and D2DGA models, respectively. The curves are collapsed onto the similarity variable $w_f=(\hat {z}/\hat {t})/\hat {w}_{0}$, representing the scaled front speed. The black and green curves in (a) represent the profiles of averaged concentration from the first and last sections of the annulus from experimental results. The green curves in (b,c) show the result of the last time step which are computed from the entire annulus. Panels (d–f) show illustrative images ($\bar {C}_{x}$) from the frontal view of the experiment in (d), the depth averaged 3D computation $\bar {c}_{x}$ in (e) and the D2DGA radially averaged concentrations $\bar {c}_r$ in (f). The times for experiment and simulations are the same. One image is taken from when the front is in the first section of the annulus and the other one from the fourth section of the annulus.

There are two additional observations from figure 5. First, the dispersion ahead of the front in figure 5(d) grows more rapidly than in the computational results. This has been observed in others of our experiments, typically at lower flow rates. Some initial mixing occurs when the gate valve is opened, but before the mixed fluid enters the annulus, which could be the cause. A similar phenomenon has been observed and explained in the work of Etrati & Frigaard (Reference Etrati and Frigaard2018), which studied miscible displacement flows in an inclined pipe with density stable configuration. Alternatively, it may be related to approaching the lower flow rate limit of the pump, where fluctuations may grow. Note that apart from controlling density, large $b$ in our experiments is controlled by setting a low flow rate.

The second observation is that in terms of predicting the changing front velocity and concentration distribution along the annulus, the D2DGA model produces more dispersion than the 3D model in this particular case. We see in figure 5(d) that the dispersive spikes are present in the annular gap and grow longer in § 4 of the annulus. However, compared with figure 4(d), they appear to extend less along the annulus and to be thinner. This feature should not be suppressed in the 3-D simulation and is still visible on close inspection. It could be that dispersive spikes are too small to be fully represented by the 3-D mesh. However, they influence the D2DGA approach directly as they are implicitly modelled. This may explain why the D2DGA model disperses more.

3.1.2. Dispersive or non-dispersive

Intuitively, an unsteady flow is always accompanied by a significant degree of dispersion; indeed unsteadiness is a form of dispersion. Many unsteady flows arise due to eccentricity (Pelipenko & Frigaard Reference Pelipenko and Frigaard2004c), but here we consider only $e=0$, so is dispersion still inevitable and also for steady front flows? According to our visual observations, in experiments and 3-D simulations, there is dispersion on the scale of the annular gap, which leads to a smearing of 2-D images of the flows. However, we have seen very clear steady fronts as well, so the question is how to identify a threshold that can be used to define different regimes and degrees of dispersion. We first present a steady front flow (figure 6) with smaller buoyancy number (EXP 48) and compare it with the case introduced previously (figure 5) to explore these differences. Ideally, we would like to quantify the dispersive effects in a way that is distinct from unsteadiness.

Figure 6.

Case EXP 48 with $(e, b, m, Re) = (0, 30, 0.8, {144})$. Panels (a–c) plot $\bar {C}$, $\bar {c}$ and $\bar {c}$ from experiment, 3-D simulation and D2DGA models, respectively. The curves are collapsed onto the similarity variable $w_f=(\hat {z}/\hat {t})/\hat {w}_{0}$, representing the scaled front speed. The black and green curves in (a) represent the profiles of averaged concentration from the first and last section of the annulus from experimental results. The green curves in (b,c) show the results of the last time step which are computed from the entire annulus. Panels (d–f) show illustrative images ($\bar {C}_{x}$) from the frontal view of the experiment in (d), the depth averaged 3-D computation $\bar {c}_{x}$ in (e) and the D2DGA radially averaged concentrations $\bar {c}_r$ in (f). The times for experiment and simulations are the same. One image is taken from when the front is in the first section of the annulus and the other one from the fourth section of the annulus.

Similar to figure 5, we find that the convergence of the experimental profiles is incomplete in the first section of the annulus, but not the last. The green lines in the computational results (figure 6b,c) again represent the last time steps calculated. Although there are significant variations in $w_f(\bar {c})$, both for low and high $\bar {c}$, the main front (at intermediate values of $\bar {c}$) evolves steadily and the criterion (3.1) is satisfied. However, we note that the main front has speed $w_f > 1$.

Second, compared with figure 5, we see significantly more dispersion. This is signalled in particular by the part of the converged $w_f$ profile at low $\bar {C}$. Rather than referring to this again as a ‘spike’, to avoid confusion, we call this feature a leading front: meaning the part of the front that advances ahead of the main front, which we threshold at $\bar {C} < 0.15$ to distinguish from the lower limit of the main front. Figure 6(c) marks this feature with a blue dashed rectangle. It corresponds to the region with light grey or yellow colour in the illustrative images in figure 6(d–f). There are always larger values for the first displacing fluids, but in cases such as in figure 5, these features are barely noticeable, occurring in the range of measurement error of $\bar {C}$ for the experiments. The leading front in figure 6 is however more substantial and signifies fluid dispersion ahead of the main front, as seen clearly in comparison to figure 5 for all three methods of approximating. Also notable is that the leading front has significantly increased in the fourth section compared with the first. On close inspection of the concentration images, we can also see that the spikes in concentration at the local gap scale (marked with circles) are significantly longer here than in figure 5.

At large $\bar {c}$, the front speed lags behind the main front speed quite significantly in the computed profiles, but less so in the experiment. This part of the curve represents the slower displacement of the wall layers around the annulus. This is calculated effectively in 3-D and is implicitly part of the fluxes in the D2DGA model. However evidently, visualizing the experiment from a side view is less effective at capturing these particular flow features. This is not surprising as in the experimental annulus, we have no means of ‘back-lighting’ the inner pipe to allow a more objective quantification of light absorption in the fluid layers between pipe and outer wall.

To have a better understanding of what is happening within the annular gap, we explore a slice along the annular gap, at fixed $\hat {x}=0$, using the 3-D results for both EXP 48 and EXP 96; figure 7. In the colour maps, $c=1$ represents the pure displacing fluid and the vectors show the normalized velocity. A long thin dispersive spike and obvious residual layers close to the walls can be observed in figure 7(a). It is the side view of such features that we see in the experimental images at the two sides of the annulus. The displacement in figure 7(b) is much more piston-like. Note that the axial velocity has a Poiseuille-type distribution that is observable both upstream and downstream in figure 7(a), away from the front. In the second case however, a more uniform velocity distribution can be observed in every axial position (figure 7b), which is due to secondary flows driven by the strong buoyancy force.

Figure 7.

Displacement flow near the front for (a) case EXP 48 $(e, b, m, Re) = (0, 30, 0.8, {144})$; (b) case EXP 96 $(e, b, m, Re) = (0, 750, 0.2, {29})$. The colour map presents concentrations across the channel gap where $\hat {x}=0$ ($c=1$ represents pure displacing fluid). The vector field shows the axial velocity.

Conservation properties of the 3-D code are generally good; see Sarmadi et al. (Reference Sarmadi, Renteria and Frigaard2021). Numerical diffusion is limited generally to 1–2 cells, e.g. see the thickness of the diffuse layer in figure 7(a), close to the wall, which is fairly characteristic. However, these thin diffuse layers are effectively dispersed in parts of the flow where the velocity difference is significant, as seen further downstream in figures 7(a) or 7(b). This is not an effect of neglecting molecular diffusion (high-Péclet-number approximation). The numerical effects would dominate physical diffusion. Taking values for water and the typical pumping speeds in the apparatus, the diffusive length scale at the exit of the annulus is $\sim$1 mm at most. As observed, the dispersive effects in the streamwise direction are significantly larger than this. However, it is possible that there could be a physical effect on the wall layer removal process. A finer mesh would better represent the spike in figure 7(a).

We have observed quite different gap-scale velocity fields associated with the different dispersion patterns. Although worthy of deeper study, the aim here is mainly to quantify different dispersion effects. We have seen in figure 6 that the pattern of dispersion ahead and behind the main front, when the displacement is steady front, has a different character. This is captured in the converged profiles $\bar {c}(w_f)$, which may equally be plotted as $w_f(\bar {c})$. On subtracting the mean pumping speed, we may define the (dimensionless) relative velocity: $w_r(\bar {c}) = w_f(\bar {c}) - 1$. We divide the relative velocity field into two components: $w_{r}^{+}$ and $w_{r}^{-}$, corresponding to the positive and negative parts of $w_{r}$, respectively. The values of $\bar {c}$ associated with $w_{r}^{+}$ characterize the dispersion downstream, moving faster than the mean flow front. Those $\bar {c}$ associated with $w_{r}^{-}$ characterize the dispersion upstream of the mean flow front, e.g. with slow removal of residual wall layers. We quantify the degree of dispersion by calculating the standard deviation at late times for each case:

(3.2)\begin{gather} \sigma_{w_{r}^{+}}=\sqrt{\frac{1}{n_{+}}\sum_{0<\bar{c}<1} (w^{+}_{r}(\bar{c}))^{2}}, \end{gather}

(3.3)\begin{gather}\sigma_{w_{r}^{-}}=\sqrt{\frac{1}{n_{-}}\sum_{0<\bar{c}<1} (w^{-}_{r}(\bar{c}))^{2}}. \end{gather}

Here, $n_+$ and $n_-$ indicate the total number of points involved in the summation calculations. Although the above calculation could be expressed as integrals, in practice, we treat the computed results in the same way as the experimental, segmenting $\bar {c}$ and summing over a finite number of discrete $\bar {c}$ values.

The results obtained are listed in table 3. As observed earlier, the side view in the experiments does not allow proper quantification of the upstream residual wall layers. Thus, only $\sigma _{w_{r}^{+}}$ is computed for the experiments, and uses data from the fourth section of the annulus ($3.6m\le \hat {z}\le 4.8m$), late in the experiment. The 3-D and D2DGA results are cleaner and clearer, so we use the entire annulus for evaluating $\sigma _{w_{r}^{+}}$ and $\sigma _{w_{r}^{-}}$. According to table 3, the dispersion along the whole annulus in EXP 48 is greater than that in EXP 96, as quantified by all methods. There are discrepancies between $\sigma _{w_{r}^{+}}$ calculated from the different models and experiment, with the experimental values the largest. However, the order of magnitude is similar and exposes differences between the two experiments. Similarly for comparisons between $\sigma _{w_{r}^{-}}$. We have also tested many other cases, not shown here for reasons of brevity. As a result, we conclude that $\sigma _{w_{r}^{+}} \le 0.08$ can be used as a threshold value to identify a non-dispersive steady front flow.

Table 3.

Standard deviation $\sigma _{w_{r}^{+}}$ calculated from experimental and computational results, and $\sigma _{w_{r}^{-}}$ calculated from 3-D and D2DGA results at a later time step, for case EXP 48 with $(e, b, m, Re) = (0, 30, 0.8, {144})$ and case EXP 96 with $(e, b, m, Re) = (0, 750, 0.2, {29})$.

However, even when $\sigma _{w_{r}^{+}} > 0.08$, the picture is not so clear. In cases such as EXP 48, the main front speed is $w_f > 1$. Thus, the computation of $\sigma _{w_{r}^{+}}$ includes the average velocity difference between the front and the mean flow, plus any effects of the dispersive spike (from values at low $\bar {c}_{y}$). Interpretation of $\sigma _{w_{r}^{-}}$ is less ambiguous, in relating only to displacement of the wall layers. To distinguish between a main front speed ($w_f > 1$) and the dispersive effect of the spike, we calculate the areas of the regions enclosed by $w_{r}^{+}(\bar {c})$ and $w_{r}^{-}(\bar {c})$, again evaluated at late times in the displacement. We define two parameters $|\bar {w}_{r}^{+}|$ and $|\bar {w}_{r}^{-}|$ as follows:

(3.4)\begin{gather} |\bar{w}_{r}^{+}| =\int_{0}^{1} w_{r}^{+} \,{\rm d} \bar{c}, \end{gather}

(3.5)\begin{gather}|\bar{w}_{r}^{-}| ={-} \int_{0}^{1} w_{r}^{-} \,{\rm d} \bar{c}. \end{gather}

These give an indication of the total dispersive area ahead of the mean position of the front, and the residual fluid left behind, respectively. We can compute these with either model, but have not found distinct differences between D2DGA and 3-D results. Whereas the standard deviations are influenced by large deviations, the above measures are not. This allows us to differentiate between effects of a steady main front ($w_f \gtrsim 1$) and dispersive spikes in the profile of $w_f$, as discussed with respect to figure 6.

To illustrate, we show a case (EXP 71) at intermediate $b=150$; see figure 8. The concentration colourmaps suggest that this case is non-dispersive. However, the spike-type profile of $w_f$ at low concentrations results in a significant $\sigma _{w_{r}^{+}}$, resulting from peak values of $w_{r}^{+}$ close to 0.13, when $|\bar {w}_{r}^{+}|$ is quite modest. In table 4, the values calculated from the D2DGA results are listed and compared with EXP 48. A first observation is that the values of $|\bar {w}_{r}^{+}|$ and $|\bar {w}_{r}^{-}|$ are similar. Second, we see that although $\sigma _{w_{r}^{+}}$ is comparable between these cases, $|\bar {w}_{r}^{+}|$ is much smaller in EXP 71. Unlike EXP 48, the main front in EXP 71 is moving at a speed much closer to the mean speed, which means that $|\bar {w}_{r}^{+}|$ is reduced (with the main contribution coming from the dispersive spike)

Figure 8.

The (a) 3-D and (b) D2DGA results of relative velocity ${w}_r=((\hat {z}/\hat {t})-\hat {w}_{0})/\hat {w}_{0}$ of case EXP 71 with $(e, b, m, Re) = (0, 150, 0.2, {145})$; (c) 3-D depth-averaged and (d) D2DGA concentration colourmap in the first and last sections of the annulus. Green lines represent the trend computed from a later time step ($\hat {t}=0.8*\hat {z}(4.8m)/\hat {w}_{0}$).

Table 4.

Areas $|\bar {w}_{r}^{+}|$ and $|\bar {w}_{r}^{-}|$ calculated from D2DGA results at a later time step for case EXP 48 with $(e, b, m, Re) = (0, 30, 0.8, {144})$ and case EXP 71 with $(e, b, m, Re) = (0, 150, 0.2, {145})$.

On exploring different cases, we have determined that $|\bar {w}_{r}^{+}|= 0.05$ is a reasonable threshold value for evaluating dispersion. To summarize, a dispersive flow is classified as such if both the following conditions are met:

(3.6a,b)\begin{equation} \sigma _{w_{r}^{+}} > 0.08 \quad \mbox{and}\quad |\bar{w}_{r}^{+}| > 0.05. \end{equation}

The flow is classified as non-dispersive otherwise. Essentially, the above quantifies the degree of smearing of the front and that it is not solely coming from the dispersive spike that is observed on the gap scale.

3.1.3. Effects of dimensionless parameters

In summary, we have now established threshold criteria to determine whether a flow is steady front, unsteady, dispersive or non-dispersive. In this section, we will classify our cases with $e=0$ and analyse further the effects of various dimensionless parameters. We also mention that some of our experimental data sets, particularly those with extremely high or low buoyancy numbers, can be classified visually. We only need to further analyse cases with moderate buoyancy force or if it is difficult to tell merely from the images. Having said this, as seen above, the values of the standard deviations and of $|\bar {w}_{r}^{\pm }|$ are also of interest for quantitative interpretation. We begin with the buoyancy number, which is often the dominant factor in primary cementing.

Keeping the viscosity ratio constant at $m=0.2$, we select cases with varying buoyancy numbers. Note that the Reynolds number often also varies with the buoyancy number, as our main control of $b$ is via the flow rate. The classification is presented in figure 9. In general, we find good agreement for the classification results obtained from all methods, except for the case with zero buoyancy. From figure 9(a), it is evident that the flow becomes more steady as the buoyancy number $b$ increases, and a magnitude of $b \sim O(10)$ is sufficient to provide a steady front in concentric annulus. It is apparently detrimental to the displacement flow if the buoyancy force is too small between fluid pairs. In addition, the discrepancy between different methods is negligible when $b > 100$, indicating that both 3-D and D2DGA models could accurately predict the steady front pattern.

Figure 9.

Classification for some of the designed cases with increasing buoyancy number, obtained from experimental and computational results as indicated by different symbols. The viscosity ratio is maintained constant at $m=0.2$. (a) Steady front or unsteady flows determined by ${\rm \Delta} w_f$; (b) non-dispersive flows determined by $\sigma _{w_{r}^{+}}$; (c) dispersive or non-dispersive flows determined by $|\bar {w}_{r}^{+}|$.

With regards to the classification of dispersive and non-dispersive flows, it is notable from figure 9(b) that the relative front velocity can be 10 % greater than the mean flow velocity $\hat {w}_0$, even when the buoyancy number is beyond 100. A piston-like displacement could be obtained only if the buoyancy number is increased to a very high level ($b \gtrapprox 10^3$). This confirms that the dispersive spike regime on the gap scale is common, but the quantity of dispersive fluids could be effectively limited if moderate buoyancy force is attained. As figure 9(c) shows, the flow is in the non-dispersive regime when $b> 100$, since the proportion of displacing fluids with faster velocity is small. Therefore, as a rough guideline for displacement flows in a concentric annulus, we transition to steady front flows with only minor dispersive effects for $b \gtrapprox 10^2$, with larger $b$ being beneficial. Both 3-D and D2DGA models are suitable at a design stage for accurate prediction. For most of our studies, the experimental results are more dispersive than the model results, but we have seen that part of this discrepancy is due to measurement error, i.e. the visualization system cannot detect all 3-D features across the gap.

The viscosity ratio $m$ appears to have a secondary role in cementing displacement flows, as described by Renteria & Frigaard (Reference Renteria and Frigaard2020), Zhang & Frigaard (Reference Zhang and Frigaard2021, Reference Zhang and Frigaard2022), Zhang et al. (Reference Zhang, Jung, Renteria and Frigaard2022). Here, we have chosen two iso-dense cases with only differing viscosity ratios ($m = 0.2,2$), the Reynolds number is maintained at the same value and $e=0$. We only focus on computational results for a clearer analysis. We know from the previous papers that the viscosity ratio mainly affects the degree of dispersion, so we plot the spatio-temporal colour map from D2DGA results in figure 10(a,b) to provide a more intuitive understanding. The region between the two dashed lines represents the concentration range from 0.1 to 0.9. In addition, the normalized front velocity is presented in figure 10(c,d). We also provide 3-D results as comparisons in the insets.

Figure 10.

The D2DGA results of (a) spatio-temporal plot of concentration of case EXP 5 with $(e, b, m, Re) = (0, 0, 0.2, 50)$; (b) spatio-temporal plot of concentration of case EXP 8 with $(e, b, m, Re) = (0, 0, 2.0, 50)$; the region within two dashed lines represents the concentration range from 0.1 to 0.9; (c) $\bar {c}$ versus $w_f$ of case EXP 5; (d) $\bar {c}$ versus $w_f$ of case EXP 8. The 3-D results of $\bar {c}$ versus $w_f$ are shown in the inset for comparison.

Focusing on the area enclosed by the dashed lines, the first case ($m=0.2$) shows a narrower region but more uniform variation of concentration. For the second case, although there is a clearer sharp front for low values of $\bar {c}_{\hat {y}}$, the wider region between the dashed lines indicates that a greater quantity of fluids are dispersing: both cases are unsteady and dispersive if we calculate ${\rm \Delta} w_f$ and $\sigma _{w_{r}^{+}}$. However, the area deviating from the mean pumping velocity is larger for the second case, with larger viscosity ratio. Specifically, $|\bar {w}_{r}^{+}|=0.1446$ in case EXP 5 ($m=0.2$) and $|\bar {w}_{r}^{+}|=0.2214$ in case EXP 8 ($m=0.8$). We find the same trend and similar values in our 3-D results. To summarize, the viscosity ratio affects the dispersive pattern to some extent and it is not beneficial to increase to $m>1$, as is intuitive. More dispersion along the whole annulus will result in a lower displacement efficiency if the pumping time is fixed. However, we confirm that viscosity ratio effects, as studied, are still secondary for large $b$. The viscosity ratio only needs to be carefully considered and designed if $b$ is fairly low (say $b <30$).

3.2.1. $b=0$

We start from the cases with zero buoyancy force to intuitively realize the effect of eccentricity. Figure 11 shows a displacement process for three cases with varying eccentricities ($e=0.2, 0.4, 0.8$, respectively), with otherwise identical parameters. As in § 3.1, figure 11(d–f) show comparative colourmaps between experimental frontal image, averaged 3-D and D2DGA results, at a single time. Figure 11(a–c) plot $\bar {c}_{\hat {y}}$ against the front velocity $w_f$, obtained from the computational results late in the displacement.

Figure 11.

Graphs of $\bar {c}$ versus $w_f$ obtained from computational results, and displacement process in the first section of the annulus presented by both experimental and computational (depth-averaged 3-D and D2DGA) results in the same time step, for the case (a,d) EXP 18 with $(e, b, m, Re) = (0.2, 0, 0.8, 500)$; (b,e) EXP 27 with $(e, b, m, Re) = (0.4, 0, 0.8, 500)$; (c,f) EXP 45 with $(e, b, m, Re) = (0.8, 0, 0.8, 500)$. Here, W and N shown in (d–f) represent the wide side and narrow side of the annulus, respectively.

Compared with the concentric flow (very similar to figure 4), we observe that with increasing eccentricity, the fluid advances progressively further on the wide side of the annulus, leaving fluid behind on the narrow side; see figure 11(d–f), taken from the first section of the annulus. We also note that the dispersion patters and gradients of the fluids are quite similar between these two experiments. There is relatively good agreement between experimental and computational results. This behaviour also corresponds to an increasing front velocity shown in figure 11(a–c). Comparisons reveal that both the D2DGA and 3-D models predict the same level of dispersion along the annulus, which is difficult to discern with the naked eye in the camera images. The shapes/slopes of the front for the different cases are consistent for all three methods. The main discrepancy found is that depth-averaged 3-D results display a more spike-like front on the wide side, whereas D2DGA results exhibit a smooth and rounded front. This difference is also noted in the maximum values of front velocity $w_f$: approximately 1.65 for D2DGA and approaching 2 for the 3-D, as shown in figure 11(c). This may be due to inertial effects, which are missing in the D2DGA model (Zhang & Frigaard Reference Zhang and Frigaard2022).

With regards to the front velocities and classification, first let us note that the data collapse nicely to define $w_f$ indicating that with eccentricity, these vertical annulus displacements remain dominantly advective. It is obvious that the conditions of ${\rm \Delta} w_f>0.1$, $\sigma _{w_{r}^{+}}>0.08$ and $|\bar {w}_{r}^{+}|>0.05$ are all met, indicating a typical unsteady dispersive flow for all three cases, as shown in figure 11(a–c). This is perhaps obvious at the outset. For these cases, we have viscosity ratio close to 1 and no buoyancy force, i.e. what mechanism can counter the tendency of the fluids to flow faster in the wider parts of the annulus?

It is interesting to consider figure 11(a–c) for the two regimes $w_f>1$ and $w_f<1$. Looking at the regime $w_f<1$, we see that the curve $\bar {c}(w_f)$ changes from convex to concave as the eccentricity increases and the area before $w_f = 1$ expands. This expanding area of $|\bar {w}_{r}^{-}|$ suggests that more of the annulus experiences residual fluids. Visually, this will likely mean that the areas of lower concentrations will extend further around the annulus azimuthally from the narrow side, as the eccentricity is increased. Since the front speeds are reduced for large $\bar {c}$, it will take significantly longer to displace all the residual fluids.

The above observation is perhaps the most damaging from the cementing perspective. Unlike other pumping operations, one cannot simply keep pumping until the residual fluid is removed. Typically, an excess volume of 20–40 % might be pumped, but there are financial, environmental and operational costs of having large excess volumes of cement slurry returning to the surface of the well. Therefore, a key objective is to minimize the area with $w_f<1$. This is done either by the use of more effective equipment to keep the casing centred in the borehole, or use of denser cement or a slower pumping flow rate to increase the buoyancy force relative to viscous forces, inducing secondary flows. The first of these is often outside of the control of the cementing company. Consequently, it is necessary to analyse the competition between buoyancy and eccentricity.

3.2.2. Competition between b and e

In this section, we explore the competition between $b$ and $e$ by considering flow types for a series of cases with differing eccentricity and buoyancy number. Figure 12 shows $\bar {c}$ versus $w_f$ obtained from 3-D results for six cases with two sets of buoyancy number ($b=30, 150$) and three sets of eccentricity ($e=0.2, 0.4, 0.8$). The insets show the frontal view from the experiments and 3-D results, as well as the narrow side view from experiments. Here, $W$ and $N$ represent the wide side and narrow side of the annulus, respectively, and the yellow blocks are the physical supports of the apparatus, reflected by the mirrors. Comparing figure 12(a, c, e) with figure 12(c, d, f), it is evident that the front dispersion has been significantly restrained by a stronger buoyant force, and much less fluid is left attached to the walls or in a narrow-side channel (for $e=0.8$). In terms of the comparisons between 3-D and experimental results, they generally display comparable front shape and position. Moreover, the 3-D results provide a clearer concentration gradient, which is hard to evaluate from greyscale experimental images.

Figure 12.

Graphs of $\bar {c}$ versus $w_f$ obtained from 3-D results, and displacement process in the first section of the annulus presented by both experimental and depth-averaged 3-D results in the same time step, for the cases (a,b) EXP 53 with $(e, b, m, Re) = (0.2, 30, 0.8, 144)$, EXP 78 with $(e, b, m, Re) = (0.2, 150, 0.8, 29)$; (c,d) EXP 58 with $(e, b, m, Re) = (0.4, 30, 0.8, 144)$, EXP 83 with $(e, b, m, Re) = (0.4, 150, 0.8, 29)$; (e,f) EXP 68 with $(e, b, m, Re) = (0.8, 30, 0.8, 144)$, EXP 93 with $(e, b, m, Re) = (0.8, 150, 0.8, 29)$. Green lines represent the trend computed from the last time step. The insets show frontal view from experimental and 3-D results, as well as the narrow side view from experiments. Here, $W$ and $N$ represent the wide side and narrow side of the annulus, respectively, yellow blocks stand for the physical supports reflected by the mirrors.

Focusing on cases with varying eccentricities but the same buoyancy number, it is difficult to discern any clear discrepancy from the changes in the front velocity, with the exception of figure 12(e), which has larger areas of both $|\bar {w}_{r}^{-}|$ and $|\bar {w}_{r}^{+}|$. The residual fluid front velocities tend to zero before $\bar {c} \to 1$, which is reflected in the narrow side view from the experiments showing an obvious residual channel, even from the beginning of the annulus. Fluids on the narrow side are difficult to displace under such high eccentricity, with relatively weak buoyant force. Apparently this issue is resolved if the buoyancy number is increased significantly, as we see in figure 12(f) for $b=150$. The front becomes flatter from the frontal view of both 3-D and experimental images, while the narrow side view reveals a much thinner residual layer. The improvement is primarily the result of a stronger secondary flow induced by larger buoyancy forces.

To have a better understanding of the secondary flow in the azimuthal direction, we present in figure 13 streamlines and velocity vectors in annular cross-sectional slices obtained from 3-D results at $\hat {z}=2.8$ m at a specific time step (slightly below the position of $\bar {c}=0.5$), for the cases shown in figure 12(e,f). The size of the yellow velocity vectors reflects the magnitude of the secondary flow, normalized with the mean velocity. We observe that the secondary flow magnitude increases with $b$ and drives fluid from the wide side towards the narrow side.

Figure 13.

Streamlines and velocity vectors of annulus cross-section obtained from 3-D results in $\hat {z}=2.8$ m at a specific time step (slightly below where $\bar {c}=0.5$), for the cases (a) EXP 68 with $(e, b, m, Re) = (0.8, 30, 0.8, {144})$ and (b) EXP 93 with $(e, b, m, Re) = (0.8, 150, 0.8, {29})$; (c) and (d) present the azimuthal velocity distribution for the corresponding cases (a to c, b to d) obtained from D2DGA results at successive time steps, and the dashed lines show $\bar {c}_r=0.5$. The colourmaps represent the normalized azimuthal velocity $\hat {v}/\hat {w}_{0}$. The size of yellow velocity vectors also stand for the magnitude of secondary flow in azimuthal direction.

The azimuthal velocity distribution for the corresponding cases, obtained from the D2DGA model at successive times are presented in figure 13(c,d). The dashed lines shows the contour $\bar {c}=0.5$. The two models predict similar maximum values for both cases. The colourmaps represent the normalized azimuthal velocity $\hat {v}/\hat {w}_{0}$. We see a characteristic switching of sign of the azimuthal velocity ahead/behind the main front, i.e. displaced fluid moves towards the wide side while displacing fluid moves to the narrow side, as also observed by Malekmohammadi et al. (Reference Malekmohammadi, Carrasco-Teja, Storey, Frigaard and Martinez2010) and Zhang & Frigaard (Reference Zhang and Frigaard2022). Again the larger $b$ results in stronger secondary flows. A simple mechanistic description of how buoyancy generates the secondary azimuthal flows is given by Maleki & Frigaard (Reference Maleki and Frigaard2018). Comparing the 2-D and 3-D results in both computations, we see an azimuthal shift, towards the narrow side, of the position of the maximum azimuthal velocity. It is interesting in figure 13(c,d) that the most severe secondary flow occurs towards the corner of the front and is terminated by the front bending into a narrow side residual channel. Looking at the residual channel, we see that the tail end of the front clearly moves upwards with time in figure 13(d), but this movement is not evident for the case with lower buoyancy number (see figure 13c). The strong secondary flow effectively assists in displacing the residual fluids stuck on the narrow side.

In summary, the buoyancy force competes with eccentricity through strong secondary flow. Under moderate buoyancy number, steady fronts and non-dispersive flow could be achieved even when $e=0.8$. Figure 14 displays the flow classification of eight cases with two different buoyancy numbers ($b=30, 150$) and increasing eccentricities ($e=0, 0.2, 0.4, 0.6, 0.8$), obtained from both D2DGA and 3-D results. Interestingly, the flow classifications do not change from those of the concentric cases. Specifically, $b=30$ (i.e. a magnitude of $b \sim O(10)$) is sufficient to provide a steady front in both concentric and eccentric annuli (figure 14a). In addition, the flow is in the non-dispersive regime when $b \approx 150$ comparing figure 14(b,c), which is identical to concentric classification. Although there are perceptible discrepancies between D2DGA and 3-D results for some of the cases, both models generally predict the same flow regimes for almost all cases covered in our experiments. There is only one data point not visible from figure 14(a), which is the value of ${\rm \Delta} w_f$ computed by the D2DGA model for the case with $b=30$ and $e=0.8$. For these parameters, we can visually still observe a steady front displacing in the plot of $\bar {c}$ versus $w_f$, but with a narrower concentration range: $0.3<\bar {c}<0.6$. This leads to a relatively large value of ${\rm \Delta} w_f$ which is not displayed as we focus in on the smaller values.

Figure 14.

Classification for some of the designed cases with two different buoyancy numbers ($b=30, 150$) and increasing eccentricities ($e=0, 0.2, 0.4, 0.6, 0.8$), obtained from both D2DGA and 3-D results as indicated by different symbols. The viscosity ratio is maintained constant at $m=0.8$. (a) Steady or unsteady flows determined by ${\rm \Delta} w_f$; (b) non-dispersive flows determined by $\sigma _{w_{r}^{+}}$; (c) dispersive or non-dispersive flows determined by $|\bar {w}_{r}^{+}|$.

Another observation is that there is relatively little sensitivity to the eccentricity for the parameters ${\rm \Delta} w_f$, $\sigma _{w_{r}^{+}}$ and $|\bar {w}_{r}^{+}|$ at the larger value of $b$. At $b=30$, we do see that $\sigma _{w_{r}^{+}}$ increases with eccentricity, as seems reasonable because we might expect that the dispersion is larger in a larger gap. When the flow is in the dispersive regime, the channelling flow gets much more severe if the annulus is more eccentric. However, at large $b$, this dispersive behaviour is suppressed; see figure 14(b) solid symbols.

We now quantify the displaced fluid that is moving slower than the mean flow velocity ($w_f<1$), which we have measured via $|\bar {w}_{r}^{-}|$, for the eight cases obtained from D2DGA and 3-D results. These are shown in figure 15. We can see an increase in $|\bar {w}_{r}^{-}|$ as the eccentricity grows, as is intuitive. The D2DGA model always predicts a larger value than the 3-D model, but within a reasonable margin of error. It appears that the classification developed for concentric annuli works well also for eccentric.

Figure 15.

Percentage of residual fluids left behind ($w_f<1$) $|\bar {w}_{r}^{-}|$ for the cases with two different buoyancy numbers ($b=30, 150$) and increasing eccentricities ($e=0, 0.2, 0.4, 0.6, 0.8$), obtained from (a) D2DGA and (b) 3-D results as indicated by different symbols. The viscosity ratio is maintained constant at $m=0.8$.

Finally, earlier we have observed distinct channels of residual fluid on the narrow side for $e=0.8$ and particularly with $b=30$. Whether we could use a threshold value of $|\bar {w}_{r}^{-}|$ to signal a narrow side channel is not clear from looking at figure 15. Here we deal with Newtonian fluids, where the narrow side channels are displaced, albeit slowly. In reality, the drilling fluids have also a yield stress and the more serious problem to consider is removing static channels of fluid on the narrow side of the annulus. Therefore, we feel that defining a threshold on $|\bar {w}_{r}^{-}|$ should wait until the non-Newtonian flows are considered.