Control of radial miscible viscous fingering

TL;DR

This paper analyzes the stability of radial viscous fingering in miscible fluids, revealing how initial radius, advection, and diffusion influence instability, and proposes a control measure based on initial radius with experimental validation.

Contribution

It introduces a stability criterion for radial viscous fingering based on initial radius, Péclet number, and log-mobility ratio, combining linear stability analysis, nonlinear simulations, and experiments.

Findings

Stability boundary approximated by M = α(r₀) Pe^{-0.55}.

Instability decreases with larger initial radius r₀.

Experimental results qualitatively agree with numerical predictions.

Abstract

We investigate the stability of radial viscous fingering (VF) in miscible fluids. We show that the instability is decided by an interplay between advection and diffusion during initial stages of flow. Using linear stability analysis and nonlinear simulations, we demonstrate that this competition is a function of the radius $r_0$ of the circular region initially occupied by the less viscous fluid in the porous medium. For each $r_0$, we further determine the stability in terms of P\'eclet number ($Pe$) and log-mobility ratio ($M$). The $Pe-M$ parameter space is divided into stable and unstable zones--the boundary between the two zones is well approximated by $M =\alpha(r_0) Pe^{-0.55}$. In the unstable zone, the instability is reduced (enhanced) with an increase (decrease) in $r_0$. Thus, a natural control measure for miscible radial VF in terms of $r_0$ is established. Finally, the results are validated by performing experiments which provide a good qualitative agreement with our numerical study. Implications for observations in oil recovery and other fingering instabilities are discussed.

Figures

Equations

$$\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}=-\frac{1}{\mu (c)}\nabla p,\mu (c)=e^{(1-c)M},$$

$${\partial }_{t}c+\mathbf{u}\cdot \nabla c=\frac{1}{Pe}{\nabla }^{2}c,$$

$$c(\boldsymbol{x},t=0)=\left\{\begin{array}[]{ll}1,&0\leq|\boldsymbol{x}|^{2}\leq r_{0}^{2}\\ 0,&\mbox{Otherwise}\\ \end{array}\right.,$$

Full text

Introduction–Hydrodynamic instabilities are ubiquitous to transport in porous media. Viscous fingering (VF) is one of these instabilities, which is observable while displacing a less mobile fluid by another more mobile fluid through porous media, and it plays critical roles in enhanced oil recovery through miscible flooding/solvent drive Lake (1989), chromatography separation Guiochon et al. (2008), pattern formation Li et al. (2009), medicines Bhaskar et al. (1992), CO2 sequestration Moortgat (2016); Amooie et al. (2017); Li et al. (2019), diffusion-limited aggregation Witten and Sander (1981), mixing Jha et al. (2011), and bacterial colonies Callan-Jones et al. (2008).

While advection is necessary for VF, diffusion stabilizes this instability in miscible systems. For a given pressure gradient, the advection velocity is uniform in a rectilinear flow, whereas, in a radial flow, the velocity is inversely proportional to the radial distance from the point of fluid injection. Effects of diffusion on stabilization of miscible fingering instabilities both in the linear and nonlinear regimes have been well understood mainly in the context of rectilinear displacement flows Homsy (1987); Pramanik and Mishra (2015). In radial VF, the effects of diffusion at the initial stages of the displacement have been studied by Tan and Homsy (1987), who concluded that the dispersion is strong enough to completely suppress the instability when the Péclet number $\sim O(1)$; otherwise the displacement is always unstable. On the other hand, Chui et al. (2015) observed shut down of the overall flow instability emerging from the dominance of diffusion over advection at the later stages of the flow. Most recently, diffusion-driven transition between two regimes of VF is captured Videbæk and Nagel (2019). Thus, there are many facets of the competition between advection and diffusion in radial flows. Contrary to this, Bischofberger et al. (2014) claim that viscosity ratio sets the velocity of the interface and three regimes of instability are obtained, for which the effects of diffusion are irrelevant. In this paper, we show that the diffusion can never be neglected when dealing with miscible fluids. Initial stable displacement is the characteristic of dominant diffusion, which is equilibrated by advection at a later stage that is identified as the transition to an unstable state dominated by advection.

Further, we ask: Can we control the competition between advection and diffusion to suppress the miscible radial VF? Diffusion, being an inherent property that depends on the displacing and displaced fluids, is difficult to tune; however, the advection can be suitably modified. Many studies focussed on controlling VF Dias et al. (2012); Zheng et al. (2015); Yuan et al. (2019) utilise time-dependent strategies to control advection. However, we achieve the same by merely modifying the initial configuration of radial displacement flow. We consider different initial finite volume of the displacing fluid in the porous medium, which is represented by an initial radial distance ($r_{0}$) of the interface from the center of the porous medium. The effects of competition between advection and diffusion on the controllability of VF are parametrised in terms of $r_{0}$ and are explained through linear stability analysis (LSA) and compared with the corresponding non-linear simulations (NLS); supported by the result of a diligently designed experiment.

Mathematical formulation and linear stability analysis–The fluids considered are Newtonian, miscible, non-reactive with $\mu_{l},\mu_{m}$ as viscosity of the less and the more viscous fluid, respectively. The non-dimensional governing equations for the flow in a two dimensional (2-D) homogeneous porous medium are constituted by the Darcy’s law and the transport equation for the solute concentration $c$,

[TABLE]

where $p$ is the hydrodynamic pressure and $\boldsymbol{u}=(u,v)$ is the Darcy velocity vector. We use $t_{f}$, the total time of fluid injection, as the characteristic time, and $\sqrt{Qt_{f}}$ as characteristic length, where $Q$ is the gap-averaged flow rate. Consequently, we obtain two non-dimensional parameters: Péclet number, $Pe=Q/D$, and the log-mobility ratio $M=\ln(\mu_{m}/\mu_{l})$, where the molecular diffusion coefficient $D$, of $c$ in the solvent fluid is assumed to be a constant.

We consider a 2-D square domain $\Omega=[-1.5,1.5]\times[-1.5,1.5]$ in the cartesian coordinates with the origin as the source of the less viscous fluid [see Appendix A.1 for computational domain]. Initial condition associated with above equations is $\boldsymbol{u}(\boldsymbol{x},t=0)=\boldsymbol{x}/(2\pi|\boldsymbol{x}|^{2})$ and

[TABLE]

where $\boldsymbol{x}=(x,y)$, and $r_{0}$ is the non-dimensional radius of initial circle occupied by less viscous fluid.

First, we perform LSA to identify the effects of diffusion at the initial stages of the displacement and the onset of VF. Assume the base state velocity to be $\boldsymbol{u}_{b}(\boldsymbol{x})=\boldsymbol{x}/(2\pi|\boldsymbol{x}|^{2})$, and the base state concentration $c_{b}$ is the solution of Eq. (3) for the initial condition (4). Analytical solution for this initial-boundary value problem is not attainable. We use the method of lines to numerically compute $c_{b}$. Spatial derivatives are discretised using sixth-order compact finite differences Lele (1992) and the resulting initial value problem is solved using a third-order Runge Kutta method. The velocity field is solved in the form of stream-functions $\psi(x,y)$, defined as $u=\partial_{y}\psi,\;v=-\partial_{x}\psi$. We introduce an infinitesimal perturbation ($\psi^{\prime},\;c^{\prime}$) such that $\psi=\psi_{b}+\psi^{\prime},c=c_{b}+c^{\prime}$ ($\left\lvert\psi^{\prime}\right\rvert,\left\lvert c^{\prime}\right\rvert\ll 1$), where $\psi_{b}$ is the stream-function corresponding to $\boldsymbol{u}_{b}$. The corresponding linearized equations are solved for $\psi^{\prime},\;c^{\prime}$ using the hybrid compact finite differences and the pseudo-spectral method. No flux boundary condition for $c^{\prime}$, and $\psi^{\prime}=0$ are used at the outflow boundary. The present LSA works as an alternative approach to study time-dependent linear system arising in miscible VF. Our interest does not lie in the wavelength selection like many other LSA Hota et al. (2015a), but to capture initial diffusion and its effect on the onset of instability. However, it must be noted that our LSA is also applicable for wavelength selection Hota et al. (2015b).

Recall that both $c_{b}$ and $c^{\prime}$ evolve temporally. Therefore, the growths of $c^{\prime},\;\psi^{\prime}$ are relative to that of $c_{b},\;\psi_{b}$, and we quantify them at each instant of time. We define the energy ratio $R(t)=E_{K}(c^{\prime},\boldsymbol{u}^{\prime})/E_{K}(c_{b},\boldsymbol{u}_{b})$, where $E_{K}(c,\boldsymbol{u})=\int\limits_{\Omega}\left[c^{2}+\boldsymbol{u}^{2}\right]d\Omega$ represents the kinetic energy at time $t$.

The energy amplification, $E(t)=R(t)/R(t=0)$ is the ratio of the energy at time $t$ to its value at $t=0$ Matar and Troian (1999). The nature of amplification as a function of time decides stability. An increasing (decreasing) $E(t)$ indicates a relative growth (decay) of the perturbations, while the presence of an extremum, if any, is of special importance. A minimum indicates transition from diffusion-dominating region to an advection-dominated one, which eventually implies the triggering of instability; while a maximum exhibits transient growths Hota et al. (2015a). Figure 1 shows the natural logarithm of energy amplification as a function of time. Evidently, $\ln(E(t))$ is a non-monotonic function with a minimum occurring for each $r_{0}>0.1$. The minimum captures the competition between advection and diffusion in the linear regime. Up to the point of minimum, the disturbances are stabilized by the diffusive base state. Larger the $r_{0}$, later the minimum is obtained, indicating a delayed onset due to the dominance of diffusive forces for a longer time. Similar qualitative results have been verified for various $M$ and $Pe$. Therefore our LSA captures controllability of VF due to the competition between the two forces.

Nonlinear simulations–Next, we perform NLS to support the estimates of LSA that the instability is delayed as $r_{0}$ increases. We solve the coupled non-linear equations, Eqs. (2)–(4), using a hybrid scheme based on compact finite differences and pseudo-spectral methods. This method has been extensively used to study instabilities in porous media (Chen et al., 2010; Sharma et al., 2019, and Refs. therein). We perform NLS for five different values of $r_{0}$, and for each $r_{0}$, we consider $M\in[0,3]$ and $Pe\in[500,10^{4}]$. A comparison of the concentration contours at a given time for the two radii in Fig. 2(a) clearly depicts that fingering instability is abated as $r_{0}$ increases. Thus, the NLS support the instability control predicted by LSA. For each simulation we compute the interfacial length Mishra et al. (2008), $I(t)=\int\limits_{\Omega}|\boldsymbol{\nabla}c|d\Omega$, which measures the temporal variation of the concentration gradient. For $M=0$, we find $I(t)=2\sqrt{\pi(t+\pi r_{0}^{2})}:=I_{0}(t)$ (say); in the absence of VF, $I(t)$ coincides with $I_{0}(t)$. We define $t=t_{on}$ as onset time of fingering in the nonlinear regime if $I(t)>I_{0}(t),\;\forall t$ satisfying $t_{on}\lesssim t<1$, and the corresponding parameter set ($r_{0},M,Pe$) is identified as unstable, otherwise stable. For each $r_{0}$, we summarize the instability in the $Pe-M$ parameter space; see Fig. 3 for $r_{0}=0.1,\;0.2$ and $0.3$. This indicates that for a fixed $M$, there is a critical Péclet number (or, similarly for a fixed $Pe$, there is a critical log-mobility ratio) for the occurrence of instability. The existence of a critical parameter for fingering instability in radial source flow with point source (i.e., $r_{0}=0$) is already identified using LSA Tan and Homsy (1987) and experiments Bischofberger et al. (2014); Videbæk and Nagel (2019). Furthermore, from NLS we observe that $t_{on}$ increases with $r_{0}$ and the stable region spans over a larger range of both $M$ and $Pe$ [see Appendix A.2 for density plots of concentration]. It is noteworthy that the boundary between the stable and unstable regions follow a scaling relation $M_{\ast}=\alpha(r_{0})Pe_{\ast}^{-\beta}$, with $0.52\lesssim\beta\lesssim 0.59$. We observe that the parameter pair ($M_{\ast},Pe_{\ast}$) lying on the boundary between the stable and the unstable regions can be well approximated by the relation $M_{\ast}=30(10r_{0}+1)Pe_{\ast}^{-0.55}$ [inset in Fig. 3]. Therefore, using this scaling relation we can approximate the stability of radial flows in homogeneous porous media.

Experiments–Our numerical results are further validated through state-of-the-art experiments. The less viscous fluid is injected at a constant flow rate $Q_{1}$ ml/s for a period of $t_{1}$ s in the Hele-Shaw cell initially filled with the more viscous fluid [see Appendix B.1, B.2 for details of experimental set up]. As soon as the required dimensional radius $\bar{r}_{0}$, is reached, the flow rate is increased to $Q^{\ast}$ and the less viscous fluid is continuously injected at this flow for a final time fixed for all the experiments. We repeat a series of experiments with $Q^{\ast}\in[0.05,0.5]$ ml/s, $M\in[0,2]$, and capture suppression of fingering instability for different values of $Q^{\ast}$ and $M$ qualitatively similar to NLS. No instability is observed for many $Q^{\ast}$ suggesting that there always exists a stable zone for each $\bar{r}_{0}$ as predicted by NLS for a range of $Pe$ and shown in Fig. 3. We use in-built MATLAB command imcontour to plot the contours at the interface. For visualisation purpose, we show the contours of one quarter the experimental image in Fig. 4(a). A delayed onset and reduced fingering corresponding to $\bar{r}_{0}=20$ mm compared to $\bar{r}_{0}=15$ mm for $Q^{\ast}=0.4$ ml/s, and $M=1.5$ are evident. This supports the fact that larger the $\bar{r}_{0}$, weaker is the advection and hence the initial competition between advection and diffusion determines the instability [see Appendix B.3 for a comparison with $\bar{r}_{0}=0$].

We use IMAGEJ Schneider et al. (2012) to quantify the experimental results. In order to avoid the noise induced by the presence of the injection pipes in the images, we utilize half of the experimental domain in our analysis. The numerical study done so far predicts stable radial displacement upto some initial time due to the dominance of diffusive forces over the advective. Hence, for the validity of the control measure, the experiments should also capture the initial radial displacement in the form of a circular displacing front. Circularity is used as the measure of the extent up to which the displacement front is circular Escala et al. (2019). We define the circularity as $C(t)=2\pi A(t)/P^{2}(t)$, where $A(t)$ and $P(t)$ correspond to the area and the perimeter of the region occupied by the less viscous fluid, respectively; so that $C(t)=1$ for a semi-circle. For $P(t)$, we consider only the length of the curved surface since the diameter does not contribute to $C(t)$. Subject to experimental errors, $C(t)$ close to $1$ but constant over the frame of time implies a circular displacing front. The maiden deviation of the circularity from the constant value indicates distortions at the front and it marks the triggering of the instability. $M=0$ works as the ideal radial source flow as no instability is observed as a result of equal viscosity of the two fluids. Hence, $C(t)$ for $M=0$ is used as the reference constant value for each $\bar{r}_{0}$. $C(t)$ for $M\neq 0$ deviates from the constant value after some initial time for each $\bar{r}_{0}$ and the time of deviation is larger for larger $\bar{r}_{0}$ (see Fig. 4(b)). In other words, the instability is triggered later for a larger $\bar{r}_{0}$ [see Appendix A.1 also for the circularity calculated from NLS]. This reassures a qualitative agreement of the experiments with LSA and NLS. The proposed control measure can be used to control VF in majority of radial source flows with miscible fluids.

Discussions and conclusion–Depending upon the application, controlling fingering instabilities is of paramount importance, such as mixing can be increased by increasing the instability; whereas for an improved oil recovery or separation process, instability should be suppressed. VF in immiscible systems are controllable by modifying the geometry Bongrand and Tsai (2018); Pihler-Puzović et al. (2018), using time dependent strategies Dias et al. (2012); Zheng et al. (2015). Though all these control mechanisms may be directly inapplicable in miscible VF Huang and Chen (2015), in this paper, we investigate a stability mechanism in miscible VF based on the principle of stabilisation in immiscible flows. Using LSA and NLS we obtain a criterion to control radial miscible VF. We utilize the competition between the advective and diffusive forces to control the miscible VF in a radial flow and show that the initial position of the miscible interface (the radius of the initial circular region containing the displacing fluid) is a control parameter. No modification in the geometry or a continuous variation in the flow rate is required unlike many other earlier studies of immiscible VF Zheng et al. (2015); Dias et al. (2012). Our theoretical results are compared with experiments that are in qualitative agreement with the numerical results. The quantitative differences between the experimental and numerical results are attributed to the gap-wise variation of the flow in the Hele-Shaw experiments compared to the 2-D numerical simulations.

Time-dependent strategies, such as changing the flow rate in a rectilinear configuration Yuan et al. (2019) or varying the gap-width of the Hele-Shaw cell in a radial configuration Chen et al. (2010), are attempted numerically in the context of miscible VF. Nevertheless, these studies lack sufficient evidences about the complete suppression of instability. Moreover, these strategies are not experimentally validated. Contrary to these, the highlights of our analysis is to successfully generate a stable radial experiments in accordance to our stability analysis and numerical simulations.

LSA predicts a delayed onset of instability with increasing $r_{0}$, while NLS predict a critical log-mobility ratio $M_{*}$ upto which there is no instability, dividing the $M-Pe$ plane into stable and unstable zones for each $r_{0}$. The log-mobility ratio and Péclet number on the boundary of the stable and unstable zones scale as $M_{*}=30(10r_{0}+1)Pe_{*}^{-0.55}$. The stable zone increases with increasing $r_{0}$. Taking $r_{0}\rightarrow 0$ we approximate the critical $Pe$ for a point source radial flow, which, for a mobility ratio $3$, is $\approx 65.8$ that is of the same order as estimated from linear stability by Tan and Homsy (1987). For fixed values of $Pe$ and $M$, it is concluded that a stable displacement for a given $r_{0}$ ensures a stable displacement for all larger $r_{0}$ despite a favorable viscosity gradient. On the other hand, an unstable displacement for a given $r_{0}$ indicates a stronger and early instability for all smaller $r_{0}$ keeping $M$ and $Pe$ fixed. This accounts to weakening of advection with an increase in the distance from the source. Experiments performed depict the validity and the applicability of the proposed control strategy.

Beside helping to understand the intrinsic properties of fundamental hydrodynamic instabilities Paterson (1981); Tan and Homsy (1987); Bischofberger et al. (2014) and pattern formation Witten and Sander (1981); Li et al. (2009), our results suggest that the controllability of miscible VF in a radial configuration could be important to predict the effectiveness of enhanced oil recovery by polymer flooding Lake (1989) and in various other similar configurations.

Acknowledgment

M.M. acknowledges the financial support from SERB, Government of India through project grant no. MTR/2017/000283. S.P. acknowledges the support of the Swedish Research Council Grant no. 638-2013-9243. C.-Y.C. is thankful to ROC (Taiwan) Ministry of Science and Technology, for financial support through Grant no. MOST 105-2221-E-009-074-MY3. V.S. acknowledges 2017 NCTU Taiwan Elite Internship Program for the financial support to visit C.-Y.C.

Appendix A Linear stability analysis and non linear simulations

A.1 Computational domain and circularity

A.2 Density plots from NLS

Appendix B Experiments

B.1 Set-up

We use a radial Hele-Shaw cell with $300\times 300\times 10$ mm3 glass plates and $b=0.5$ mm gap width. A T-junction made up of a hypodermic syringe needle ($1.2$ mm diameter and $40$ mm length) bent in an L-shape and carefully embedded into a pipe, is used for filling the two fluids in the Hele-Shaw cell. The t-junction avoids the hassle of (a) changing the pipes for different fluids, and (b) having a hole in each glass plate. A syringe pump (Cole-Parmer-D201253) is used for injecting the less viscous fluid. The dynamics are captured with Sony FDR-AX40 camera.

B.2 Initial condition

The less viscous fluid is injected at a constant flow rate $Q_{1}$ ml/s for a period of $t_{1}$ s in the Hele-Shaw cell initially filled with the more viscous fluid. The initial volume of the injected less viscous fluid, $Q_{1}t_{1}$ ml, is so chosen that it leads to a stable displacement of the more viscous fluid until the invading fluid occupies a circular region of radius $\bar{r}_{0}\;$ mm. Recall that during this stable displacement the interface between the two fluids experiences diffusive spreading proportional to a length $\sqrt{Dt_{1}}$, which also contributes in $\bar{r}_{0}$; we measure $\bar{r}_{0}=\sqrt{Dt_{1}}+\sqrt{Q_{1}t_{1}/(\pi b)}$. Here, $D=10^{-9}$ m2/s is the molecular diffusion coefficient of glycerin in water D’Errico et al. (2004).

B.3 Effect of increasing $\bar{r}_{0}$ on VF