Effective Rheology of Two-Phase Flow in Three-Dimensional Porous Media: Experiment and Simulation

Abstract

We present an experimental and numerical study of immiscible two-phase flow of Newtonian fluids in three-dimensional (3D) porous media to find the relationship between the volumetric flow rate (Q) and the total pressure difference ([Formula: see text]) in the steady state. We show that in the regime where capillary forces compete with the viscous forces, the distribution of capillary barriers at the interfaces effectively creates a yield threshold ([Formula: see text]), making the fluids reminiscent of a Bingham viscoplastic fluid in the porous medium. In this regime, Q depends quadratically on an excess pressure drop ([Formula: see text]). While increasing the flow rate, there is a transition, beyond which the overall flow is Newtonian and the relationship is linear. In our experiments, we build a model porous medium using a column of glass beads transporting two fluids, deionized water and air. For the numerical study, reconstructed 3D pore networks from real core samples are considered and the transport of wetting and non-wetting fluids through the network is modeled by tracking the fluid interfaces with time. We find agreement between our numerical and experimental results. Our results match with the mean-field results reported earlier.

Keywords: Dynamical pore network model; Reconstructed porous media; Steady-state two-phase flow; Two-phase flow experiment.

Fig. 1

The experimental setup of two-phase flow of deionized water as the wetting fluid and air as the non-wetting fluid through the column of borosilicate glass beads

Fig. 2

Three sample pressure plots of the time evolution of the experiment are shown to demonstrate the chaotic and transient nature of the two-phase flow. Plots from the high-Ca regime (a), transition point (b), and low-Ca regime (c) illustrate the characteristic pressure behavior observed

Fig. 3

Plot of dimensionless excess pressure drop ($B-B_t$) as a function of the capillary number (Ca) for the nonlinear and linear flow regimes obtained from the experiment. Each data point represents the averaged steady-state pressure gradient for the respective Ca. The error bars are obtained from the fluctuations in the pressure drops in the steady state. The two scaling exponents for the low and the high-Ca regimes obtained from the slopes are $0.46\pm 0.05$ and $0.99\pm 0.02$ respectively

Fig. 4

The schematic of one link between the two nodes i and j is shown in (a). The pore space of each link is divided into three pore parts, two pore bodies at two ends and one pore throat in between. The total length ($l_{ij}$) of the link is equal to ${\mathit{\Lambda }}_1+{\mathit{\Lambda }}_2+{\mathit{\Lambda }}_3$, the sum of each pore part. The presence of multiple interfaces between the wetting (white) and non-wetting (gray) fluids in a link is shown in the bottom of (b), and the variation of capillary pressure $p_c(x)$ as a function of the interface position (x) for each interface is shown in the top of (b). $p_c=0$ at the two ends of the link and is maximum, equal to the threshold pressure $p_t$, at the middle of the tube. $p_t=4\mathit{\gamma }cos\mathit{\theta }/r_t$ is the minimum pressure required for the non-wetting fluid to invade the pore

Fig. 5

The time evolution of two-phase flow through the reconstructed networks where the wetting and the non-wetting fluids are colored by blue and red, respectively. The few links in black are the dead ends which are connected only at one node and therefore removed from the network. The three images from left to right in each row, respectively, correspond to the initial condition of the system and after 0.1 and 0.3 pore volumes of fluids have passed. The three rows from top to bottom correspond to the samples A (berea), B (sandpack) and C (sandstone), respectively. The overall flow is in the positive x direction. The periodic boundary condition is implemented in the same direction, by making a mirror image of the original reconstructed network and then connected together with the original. These two parts are shown by the two cuboids in the figures. Here, the system is initialized by filling the links with non-wetting fluid from $x=0$ until the required saturation is obtained and then filling the rest with the wetting fluid. In these figures, the non-wetting saturations are 0.3 for sample A and 0.5 for B and C

Fig. 6

Plots of non-wetting fractional flow ($F_{\mathrm{nw}}$) in the steady state as a function of non-wetting saturation ($S_{\mathrm{nw}}$) for the three different networks. The dashed diagonal straight lines correspond to $F_{\mathrm{nw}}=S_{\mathrm{nw}}$, a system of miscible fluids will follow that line. The results for two different capillary numbers, $\text{Ca}={10}^{-1}$ and ${10}^{-2}$, are shown. Notice that the curves approach the diagonal straight line for the higher value of Ca. In the inset for sample A, we plot the fractional flow as a function of pore volumes ($N_{\mathrm{v}}$) of fluids passed, where $F_{\mathrm{nw}}$ fluctuates around an average value in the steady state

Fig. 7

Total pressure drop ($\mathrm{\Delta }P$) in the steady state for the three networks as a function of $S_{\mathrm{nw}}$. The capillary number $\text{Ca}={10}^{-2}$. $\mathrm{\Delta }P$ reaches a maximum at an intermediate saturation, which is due to the increasing number of interfaces causing higher capillary barriers. In the inset of sample A, $\mathrm{\Delta }P$ is plotted as a function of the pore volumes of fluids passed, which shows the evolution of steady state

Fig. 8

Plot of overall pressure drop as a function of the capillary number in the steady state for different networks. For network A, results are shown for $S_{\mathrm{nw}}=0.3$ with two different viscosity ratios $M=1.0$ and 0.1. For B and C, two different saturations $S_{\mathrm{nw}}=0.3$ and 0.5 with $M=1.0$ are considered. Measurements of the threshold pressures for each simulation are shown in the insets of each figure, where $P_t$ is obtained from the y-axis intercepts of $\mathrm{\Delta }P$ versus $\sqrt{\text{Ca}}$ plots. Using the values of $P_t$, the scaling exponents are then obtained from the slopes of log–log plots of ($\mathrm{\Delta }P-P_t$) versus Ca