Inertia Effects in the Dynamics of Viscous Fingering of Miscible Fluids in Porous Media: Circular Hele-Shaw Cell Configuration

Abstract

We present a numerical study of viscous fingering occurring during the displacement of a high viscosity fluid by low viscosity fluid in a circular Hele-Shaw cell. This study assumes that the fluids are miscible and considers the effects of inertial forces on fingering morphology, mixing, and displacement efficiency. This study shows that inertia has stabilizing effects on the fingering instability and improves the displacement efficiency at a high log-mobility-viscosity ratio between displacing and displaced fluids. Under certain conditions, inertia slightly reduces the finger-split phenomenon and the mixing between the two fluids.

1. Introduction

Enhanced Oil Recovery (EOR) is one of the first responses to the increasing global demand for energy. However, EOR operations including water flooding and gas injection suffer from poor sweep efficiency, particularly in highly heterogeneous reservoirs. The sweep efficiency becomes even worse in the presence of heavy oils. The poor sweeping efficiency in EOR operations is the result of the development of so-called viscous fingering [1]. Understanding the fingering phenomenon is of primary importance.

The viscous fingering phenomenon goes back to the pioneering work of Saffman and Taylor in 1958 [2]. This instability occurs at the interface between two fluids when a less viscous fluid displaces a more viscous fluid and appears as finger-like patterns. Fingering instability occurs at various scales ranging from micro-scale to several kilometers in fractured reservoirs [3], in chromatography separation [4,5], Enhanced Oil Recovery [6,7], CO₂ sequestration [8,9], and pollutant removal processes [10].

The pioneering work of Saffman and Taylor is the starting point for many theoretical, experimental, and computational studies [11,12,13]. These studies explored the onset of instability or finger-like patterns, their shape and tip-splitting [14,15,16]. Studies also investigated the effects of the fluid rheology of geometrical disturbances on the development of fingering instability [17,18,19,20,21,22].

The Hele-Shaw cell is a canonical setup used to comprehend fingering instability occurring during fluid displacement in porous media. Most numerical studies were conducted in semi-finite rectilinear configurations using Darcy’s Law. These studies showed that viscous fingering was found in miscible fluid flows because of the change in effective viscosity during mixing [23,24]. The viscosity contrast between the two fluids together with a non-homogeneous spatial distribution of permeability and fluid speed control the morphology of the finger patterns [25]. A Positive log-mobility–viscosity ratio is found to produce instabilities at the rear section of the fluid interface, while a negative log-mobility–viscosity ratio generates forward-directed fingers [26]. Studies have explored the possibility of controlling fingering dynamics and morphology using a circular Hele-Shaw cell [27,28,29]. For instance, some studies explored the control of fingering instability with the variation in time of the gap between the plates constituting the Hele-Shaw cell [27,28]. Other studies examined the suppression of fingering instability using a tapered Hele-Shaw cell [29]. It was also shown that fingering instability can be controlled by replacing the upper plate in the Hele-Shaw cell with an elastic membrane. It has been found that the resulting fluid–structure interaction changes the interfacial patterns and delays the onset of fingering instabilities [30].

The effects of inertia on finger instability received little attention. Almost all of the numerical and theoretical studies were based on Darcy’s law, which dates back to 1856 [31]. Darcy’s model neglects inertial forces; it only involves viscous effects. However, inertia can become non-negligible with the increase in the flow speed. This situation can be encountered during gag injection in the reservoir. It has been found that the finger width is the result of a balance between viscous and inertial forces and that the finger width can increase with increasing inertia or flow speed [32]. To rationalize their observations, the authors in [32] conducted numerical simulations using a modification of Darcy’s law which includes inertia. The authors only found qualitative agreement with the experimental data [32]. In the literature, one can find several empirical attempts to include inertial forces in the flow in porous media [33,34]. For example, to consider inertial effects in porous media, Forchheimer included an empirical parameter β in Darcy’s equation [35]. However, finding a correct value of β to distinguish Darcy and post-Darcy regions where inertial effects are important is not a trivial task. The researcher used the Reynolds number to distinguish between these two regions [36,37,38,39,40]. Recently, a criterion to consider inertial effects was defined in terms of a Forchheimer number; $Fo\ge 0.11$. β in the definition of a Forchheimer number, which was earlier considered as independent of the flow speed, is now considered varying with the speed of the flow. This makes the application of this empirical criterion difficult [38]. Two other empirical formulations, such as the Hagen–Poiseuille and the Ergun models, were introduced. Both models rely on one or more empirical constants.

There are very few theoretical models taking into account inertial effects in the problem of fluid displacement in Hele-Shaw cells. Rabaud et al. introduced a new equation where viscous effects are added to the Euler equation and inertia is included as a drag term [41]. An inertial correction to the Darcy law in a Hele-Shaw cell was proposed by Ruyer-Quil [42], who used a perturbation method and a polynomial approximation for the velocity field. The obtained equation is optimal in the sense that every method of weighted residuals will converge to it. Using this model, Ruyer-Quil showed good agreement with the experimental data, given in Rabaud et al. [41]. To the best knowledge of the authors, there are no numerical studies that considered the effect of inertia on displacement efficiency, and fluid mixing [42,43,44,45]. Our group has been involved in multiphase flow in porous media, both at high and low interfacial velocities [46,47,48].

In this work, we investigate the dynamics of the fingering phenomenon in the case of miscible fluids within a circular Hele-Shaw cell. This work parallels, to some extent, our previous work on viscous fingering in the case of immiscible fluids [49]. The present study includes the effects of inertia on finger morphology, displacement efficiency, and fluid mixing. We used the Finite Element Method (FEM) to solve the problem at hand. To depict the effects of inertia, we used both Darcy’s and modified Darcy’s laws based on the Ruyer-Quil model.

2. Mathematical Modeling

We consider a 2D flow between two circular plates that form a circular Hele-Shaw cell, where a less viscous fluid displaces more viscous fluid. The outer cell radius is assumed to be fifty times the inlet radius of the central injection; see Figure 1. We believe that a ratio of fifty is a reasonable one. With a smaller ratio, fingers may not have enough distance to develop and reach their stable pattern. On the other hand, a large ratio will result in alengthy computational time.

Fluid I enters from the center and spreads radially to displace Fluid II towards the peripheral outlet of the flow domain or Hele-Shaw cell. We assume that both fluids are Newtonian, miscible, non-reactive, and incompressible. We also assume that the dynamic viscosity of Fluid II is a function of solute concentration c, contained in Fluid I. The solute concentration disperses in Fluid II such that $c=0$ and $c=1$ correspond to Fluid I and Fluid II, respectively. The flow configuration is assumed to be horizontal, hence the effect of gravity is neglected. The governing equations are given as follows:

$$\overrightarrow{\nabla }·\overrightarrow{u}=0$$

(1)

$$\overrightarrow{\nabla }p=-\frac{\mu \left(c\right)}{k}\overrightarrow{u}$$

(2)

$$\varphi \frac{\partial c}{\partial t}+\left(\overrightarrow{u}·\overrightarrow{\nabla }\right)c-\overrightarrow{\nabla }·\left(D_m\overrightarrow{\nabla }c\right)=0$$

(3)

Equation (1) indicates the incompressibility of the fluids. Equation (2) is Darcy’s law, where $p,\mu ,k$ indicates the pressure, dynamic viscosity, and absolute permeability of the medium, respectively. Dynamic viscosity is considered as a function of solute concentration such that [50]:

$$\mu \left(c\right)={\mu }_1{e}^{Rc}$$

(4)

Here, R is the log-viscosity ratio of the dynamic viscosities,

$$R=ln\left(\frac{{\mu }_2}{{\mu }_1}\right)$$

(5)

Equation (3) is the transport equation of solute concentration c, where ϕ indicates the porosity of the medium and D is the molecular diffusion coefficient. The concentration of the species is defined as $c=c_1={\rho }_1{S}_1$, where ${\rho }_1$ and $S_1$ are the density and fluid saturation of Fluid I, respectively. Since the effects of gravity are not included in our model, arbitrary values can be assigned for density without affecting the solution. For the sake of simplicity, we assume ${\rho }_1={\rho }_2=1$. The governing equations become,

$$\varphi \frac{\partial S_i}{\partial t}+\overrightarrow{\nabla }·\left(S_i\overrightarrow{u}\right)-D\frac{{\partial }^2{S}_i}{\partial x^2}=0$$

(6)

$$\overrightarrow{\nabla }p=-\frac{\mu \left(S_i\right)}{k}\overrightarrow{u}$$

(7)

$${\displaystyle \sum }_{i=1,2}{S}_i=1$$

(8)

where i = 1, 2 denotes Fluid I and Fluid II, respectively. Equations (6) and (7) are new forms of Equations (3) and (4), respectively. Equation (8) is the total mass conservation in terms of saturation of Fluids I and II. To include inertia, we used a modified Darcy’s law or Ruyer’s model [40]. In this case, Equation (7) is replaced by the following equation:

$$\frac{6}{5}\frac{\partial \overrightarrow{u}}{\partial t}+\frac{54}{35}\left(\overrightarrow{u}·\nabla \overrightarrow{u}\right)=-\frac{1}{\rho }\nabla p-\frac{\mu \left(S_i\right)}{k}\overrightarrow{u}$$

(9)

Equation (9) is derived using perturbation and weighted residuals methods with a polynomial approximation of the velocity field. This method is similar to the boundary-layer integral method, where a parabolic velocity profile is used to determine the boundary layer thickness. The difference between the two methods lies in the choice of test or (weighted) functions. Equation (9) is optimal in the sense that every method of weighted residuals will converge to it as the number of test functions is increased.

We use FEM-based commercial software, COMSOL Multiphysics 5.3, to solve the above systems of equations. We solved Darcy’s model with Equations (6)–(8), and modified Darcy’s with Equations (6), (8) and (9). The systems of partial differential equations were solved with the appropriate initial and boundary conditions. Initially, we assume that the whole domain is filled with motionless Fluid II. Fluid I enters from the inlet and radially spreads in the domain and displaces Fluid II. The velocity and saturation at the inlet ($r=r_{inner}$) are:

$$\begin{gathered} (\mathrm{i})\overrightarrow{n}·\overrightarrow{u}=U_{inlet} \\ (\mathrm{ii})S_2=1, \end{gathered}$$

where $\overrightarrow{n}$ is the radial normal unit vector.

At the outlet ($r=r_{outer}$), we impose constant pressure and no-flux of the saturation such that:

$$\begin{gathered} (\mathrm{iii})p=0 \\ (\mathrm{iv})\overrightarrow{n}·\nabla S_i=0. \end{gathered}$$

3. Results

In this section, we summarize the results of the simulation conducted both with Darcy’s law and modified Darcy’s law. The simulations were conducted in a COMSOL environment. The numerical model and mesh size are validated qualitatively using experimental observations by Chen et al. [51]. Figure 2 and Figure 3 summarize the mesh sensitivity. The mesh was gradually refined until the pattern became independent of the mesh size [6,52]. Figure 2 summarizes the mesh sensitivity. All the simulations were conducted using the finest mesh that allowed us to recover the pattern observed experimentally in [51]; see Figure 2.

In the simulations, we set $k=1\times {10}^{-9} {\mathrm{m}}^2$, $\varphi =0.5$ and the initial viscosity of Fluid II in the domain as ${\mu }_2=1\times {10}^{-3} \mathrm{P}\mathrm{a}.\mathrm{s}$. We performed numerical simulations with and without inertial effects after setting the inlet velocities as $U_{inlet}=2\times {10}^{-4} \mathrm{m}{\mathrm{s}}^{-1}$, $U_{inlet}=2\times {10}^{-3} \mathrm{m}{\mathrm{s}}^{-1}$ and $U_{inlet}=2\times {10}^{-2} \mathrm{m}{\mathrm{s}}^{-1}$. Similarly, to study the effects of viscosity, we performed simulations for $R=2$, $R=4$, $R=6$ and $R=8$ for all the above velocities with ${\mu }_1<{\mu }_2$. So, in total, we performed twenty-four simulations with twelve combinations of viscosity and three velocities each for Darcy and modified Darcy models. In these conditions, the value of Reynolds $\left(Re=\frac{\rho Uh}{\mu }\right)$ and Peclet $\left(Pe=\frac{Uh}{D}\right)$ numbers range from 2.2 × 10⁻⁵–2.2 × 10⁻³ and 0.015–1.5, respectively. In the following sections, we discuss the effects of inertia on fingers morphology, displacement efficiency, and mixing area.

3.1. Fingers Morphology

Figure 3, Figure 4 and Figure 5 depict finger patterns obtained for various values of Reynold, Peclet, and the log-mobility–viscosity ratio. These figures confirm that the fingering instability is related to the ratio of the viscosities of the two fluids. For instance, when the log-mobility–viscosity ratio R is equal to 2, almost no instability is found. The instability starts to develop with the increase in R. Increasing the log-mobility–viscosity-ratio results in an increase in the number of fingers and their split. However, at the onset of the stability when R is equal to 4 and there are relatively low Reynolds and Peclet numbers, inertia reduces the splitting of the fingers; see Figure 3c,d. These two figures compare the finger patterns obtained with Darcy’s law and modified Darcy’s law when R = 4, $Pe=1.45\times {10}^{-2} Re=2.2\times {10}^{-5}$ at the time $t=4000 \mathrm{s}$. Figure 3, Figure 4 and Figure 5 show that with increasing R flow the effect of inertia diminishes until it has almost no significant effect on the splitting of the fingers. The nonlinearity of the inertia term seems to saturate the growth of the fingering instability.

3.2. Displacement Efficiency

The displacement efficiency is an important factor in fluid displacement problems, particularly in Enhanced Oil Recovery (EOR) and soil decontamination. Fingering instability is responsible for the drop of sweeping or displacement efficiency during oil recovery operations. The displacement efficiency is defined as the saturation of Fluid I displacing Fluid II in the flow domain. The values of the displacement efficiency obtained for various combinations of R, Pe, and Re with Darcy and modified Darcy models are compared in Figure 6. This figure shows that the fluid displacement efficiency increases rapidly to achieve its asymptotic value (100%) for R = 2 and high values of Reynolds and Peclet numbers. This result is expected because when R = 2, no instability is found and Fluid I advances as a circular front, entirely expelling Fluid II from the Hele-Shaw cell. Increasing the log–viscosity ratio results in a decrease in the displacement efficiency because of the development of fingering instability. Figure 6 also shows that the displacement efficiency obtained with modified Darcy’s law starts to depart from the ones obtained with Darcy’s law when R increases. At R = 8, the simulation with inertia indicates a high displacement efficiency. This means that inertia reduces the effects of fingering instability. At the onset of the instability, R = 4, inertia enhances fingering instability and reduces the displacement efficiency. Inertia affecting the displacement efficiency is an important result for EOR using gas injection. This means that controlling the speed of gas injection in the oil reservoirs can improve the sweeping efficiency of the gas injection by reducing the effect of fingering instability through inertia. Some gases such as CO₂ are miscible with oil and their log-mobility–viscosity ratio, R, can be monitored in such a way that inertia alleviates fingering instability and increases the sweeping efficiency.

3.3. Mixing Area

The mixing area is another important parameter; it is of fundamental importance in miscible fluids. This parameter estimates the diffusion of Fluid I within Fluid II. The mixing area measures the ratio between the area of the mixed flow and the area of the whole flow domain. In other words, the mixing area is the total surface where the saturation is between 0 and 1. Saturation equal to 1 means only displaced fluid occupies the whole domain, while saturation equal to 0 means only displacing fluid occupies the whole domain.

A high mixing area indicates good mixing and implies strong instabilities. Figure 7 shows that the mixing area increases with time until it reaches its maximum value when fingers reach the outlet of the flow domain. At this moment, the width of the fingers starts to increase and facilitates fluid drainage through the fingers. The mixing area starts to decrease from the time at which the breakthrough of Fluid I occurs. Similarly, the displacement efficiency increases rapidly until Fluid I reaches the peripheral outlet. From this moment on, the displacement efficiency increases slowly towards its asymptotic value. Inertia does not have a significant effect on the mixing area. Figure 7 compares the mixing areas obtained at different values of R for fixed values of Reynolds and Peclet numbers. Figure 8 indicates that inertia has very little influence on the mixing areas.

4. Conclusions

We investigated the effects of inertia on the onset and spreading of fingering instability in a circular Hele-Shaw cell. We have considered miscible fluids and discussed the effects of inertia on finger morphology, displacement efficiency, and mixing area for various log-mobility–viscosity ratios, Reynolds, and Peclet numbers. Inertia is included using modified Darcy’s model proposed by Ruyer-Quil [40]. Our simulations indicate that the onset of fingering instability occurs when the log-mobility–viscosity ratio, R, is greater than 2. The simulations show that inertia can alleviate the splitting of the fingers at the onset of instability when Reynolds and Peclet numbers are small. Simulations also show that inertia improves the displacement efficiency when the log-mobility–viscosity ratio is high. However, at the onset of the fingering instability at small values of the log-mobility–viscosity ratio, inertia reduces the displacement efficiency. Hence, inertia enhances fingering instability at small values of R, while it reduces it at large values of R. Simulations indicate that inertia does not affect the mixing area. Inertia has a significant effect on the displacement efficiency when Reynold and Peclet numbers and the log–viscosity ratio are high. This is an interesting result for the simulation of gas injection EOR operations. Simulation with Darcy’s model can underestimate the recovery factor of EOR operation when log-mobility–viscosity ratios are high or overestimate it when log-mobility–viscosity ratios are small.