Geofluids

Volume 2018, Article ID 8269645, 13 pages

https://doi.org/10.1155/2018/8269645

## Fluid Interfaces during Viscous-Dominated Primary Drainage in 2D Micromodels Using Pore-Scale SPH Simulations

Institute of Applied Mechanics, University of Stuttgart, Pfaffenwaldring 7, 70 569 Stuttgart, Germany

Correspondence should be addressed to Holger Steeb; ed.tragttuts-inu.uabhcem@beets.regloh

Received 27 December 2017; Accepted 8 April 2018; Published 14 June 2018

Academic Editor: Emanuele Romano

Copyright © 2018 Rakulan Sivanesapillai and Holger Steeb. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We perform pore-scale resolved direct numerical simulations of immiscible two-phase flow in porous media to study the evolution of fluid interfaces. Using a Smoothed-Particle Hydrodynamics approach, we simulate saturation-controlled primary drainage in heterogeneous, partially wettable 2D porous microstructures. While imaging the evolution of fluid interfaces near capillary equilibrium becomes more feasible as fast X-ray tomography techniques mature, imaging methods with suitable temporal resolution for viscous-dominated flow have only recently emerged. In this work, we study viscous fingering and stable displacement processes. During viscous fingering, pore-scale flow fields are reminiscent of Bretherton annular flow, that is, the less viscous phase percolates through the core of a pore-throat forming a hydrodynamic wetting film. Even in simple microstructures wetting films have major impact on the evolution of fluid interfacial area and are observed to give rise to nonnegligible interfacial viscous coupling. Although macroscopically appearing flat, saturation fronts during stable displacement extend over the length of the capillary dispersion zone. While far from the dispersion zone fluid permeation obeys Darcy’s law, the interplay of viscous and capillary forces is found to render fluid flow within complex. Here we show that the characteristic length scale of the capillary dispersion zone increases with the heterogeneity of the microstructure.

#### 1. Introduction

Assessing the stability and evolution of saturation fronts, or, from a pore-scale point of view, interfaces between immiscible bulk fluid phases, is key with respect to understanding a multitude of subsurface processes. Examples include sequestration of carbon dioxide in geological media, groundwater contamination remediation, or enhanced oil recovery. Depending on the governing capillary number (Ca), viscosity ratio (), properties of the porous microstructure, and boundary conditions, primary drainage results in flow regimes as diverse as viscous fingering, stable displacement, or capillary finger branching [1]. In traditional macroscopic continuum models, a phenomenological extension of Darcy’s law is assumed to be applicable with relative permeability and capillary pressure functions representing constitutive model inputs [2]. Calibration of these constitutive relations in light of a particular flow regime typically renders them nonlinear and hysteretic [3–5]. In an attempt to face the latter, contemporary models acknowledge the role of interfacial areas in hysteresis [6–10] or explicitly account for mass-exchange between hydraulically reservoir-connected and reservoir-disconnected subphases [11, 12].

Considerable effort has been devoted to studying two-phase flow patterns at the length scale of pore-networks, both experimentally [13–18] and numerically [19–25], providing a reliable set of data for the - phase diagram of drainage displacement patterns as introduced by [1]. However, the pore-scale dynamics of fluid interfaces and the mechanisms by means of which they evolve remain poorly understood.

To this end, recent advances in pore-scale imaging-based characterization methods (see review [26]) that enable the fast visualization of two-phase flow at pore-scale resolution, most notably microscopy imaging of thin micromodels [18, 27–29], X-ray computed tomography [30–34], and confocal microscopy [35, 36], have provided valuable insights into the interplay of viscous, capillary, gravitational, and inertial forces constituting the complexity of interface dynamics at the pore-scale. For instance, free-energy driven Haines jumps have been confirmed as dominant displacement mechanism for flow at small capillary numbers [29, 32]. Clearly, these observations deviate from the assumptions underlying generalized Darcy flow. Besides experimental approaches, we consider direct numerical simulations (DNS) to be an important complementary tool for quantitative characterization of multiphase flow in porous media [26, 37].

Our method of choice is a quasi-incompressible Smoothed-Particle Hydrodynamics (SPH) model [38–41] which incorporates the Navier-Stokes equations together with the Continuum Surface Stress method [42, 43] to account for the interfacial balance equations. Advantages of SPH in the context of two-phase flow in porous media include its mesh-free nature. The reproducing kernel interpolation of SPH does not require nodal integration points (particles) to be distributed on grids or meshes. The latter renders spatial discretization of complex pore spaces less computationally expensive as compared to traditional grid or mesh-based methods. Typical SPH methods use an updated Lagrangian approach; that is, particles are advected in space according to their local advection velocity. Hence, nonlinear convection terms are not required to be modeled. The latter further implies that phase indicator fields are advected through particle motion. As a result, implementation of the Continuum Surface Stress method does not require interface-tracking such that SPH constitutes an attractive approach to model formation, coalescence, and fragmentation of fluid interfaces in complex pore spaces. Moreover, the updated Lagrangian formulation simplifies modeling of locally large Reynolds numbers [44, 45]. Disadvantages of SPH include its high computational costs associated with repetitive use of neighbor searching algorithms. In the context of using explicit time integration schemes, time stepping stability criteria such as the CFL-condition may be rather restrictive. Moreover, since uncorrected SPH interpolation stencils lack the Kronecker delta property, the application of essential boundary conditions remains nontrivial [46, 47].

In this work, we present pore-scale resolved DNS of two-phase flow in 2D partially wettable porous media of particulate microstructure. We perform numerical experiments of saturation-controlled primary drainage for various magnitudes of and . While much effort has been devoted to studying the pore-scale distribution of fluid interfaces near capillary equilibrium, pore-scale numerical studies for viscous-dominated flow remain comparatively sparse. Rather than studying the emerging macroscopic displacement patterns we discuss pore-scale flow fields. For viscous fingering we discuss the formation of hydrodynamic wetting films and their impact on the evolution of specific interfacial area and interfacial viscous coupling effects. For stable displacement we discuss the microscopic dispersion of stable saturation fronts (capillary dispersion zone) and the role of microstructure heterogeneity. In an attempt to elucidate the trade-off between modeling assumptions, topological constraints of 2D simulations, and physical meaningfulness, we critically discuss our numerical approach.

#### 2. Methods

##### 2.1. Direct Numerical Simulation Approach

We model isothermal flow of a Newtonian, quasi-incompressible wetting () and nonwetting () fluid phase through the pore-space () of a rigid solid matrix (). Bulk fluid balance equations for mass density and linear momentum in phase domains readrespectively. The Newtonian viscous extra stress tensor reads and denotes local velocity. Dynamic viscosity is considered constant. A superimposed dot denotes the material time derivative. Fluid pressure is related to bulk mass density by a stiff barotropic equation of state such that density fluctuations are small implying quasi-incompressibility. A constant, positive backpressure is included to avoid negative pressures that otherwise lead to a numerical tensile instability which SPH methods are prone to [48].

Immiscibility and solid wettability are accounted for using interfacial balance equations for all Gibbs material interfaces . Interfacial balance equations are expressed in terms of jump conditions using the Rankine-Hugoniot jump operator . The interfacial mass balanceswhere denotes interfacial velocity, imply kinematic coupling along the interfacial normal vector . Interfacial balances of linear momentum readwhere denotes the surface divergence operator and the interfacial Cauchy stress tensor . Interfacial tensions are considered constant which is considered reasonable in the absence of surface-active agents. Solid wettability properties are prescribed by setting parameters , , and appropriately.

An alternative expression of the interfacial balance of linear momentum (4) for the fluid-fluid interface readswhere is twice the mean curvature. The first term in (5) implies interfacial viscous coupling whereas the second term introduces a pressure jump condition according to the Young-Laplace equation [49].

Our DNS approach [41, 45, 50] is based on the Continuum Surface Stress method [42, 43] whereby bulk balance equations are formulated for the whole-fluid domain . The latter is achieved by immersing jump conditions (4) using interface Dirac delta distributions which are compactly supported on points of only. The resulting whole-fluid balance of linear momentum is expressed aswhere the solid surface identity tensor is introduced to remove the otherwise unbalanced contact line stress normal to the solid surface [45, 50], that is, to reduce contact line forces to planes tangent to the solid surface. Discretization of (6) using SPH leads to the intuitive motion equationfor a fluid particle with lumped mass . Upon discretization, internal bulk forces due to viscous diffusion and pressure gradient result in discrete interaction forces ( and ) that act between and its neighboring particles .

##### 2.2. Simulation Setup and Microstructures

The total computation domain (Figure 1) is comprised of wetting phase (WP) and nonwetting phase (NWP) reservoirs denoted by and , respectively, as well as the porous computational unit cell . The unit cell has side lengths mm and mm in direction of the unit vectors and , respectively. The microstructure exhibits -periodicity in direction of such that periodic boundary conditions are applied with respect to . The latter avoids boundary artifacts that otherwise arise due to wall confinement. Initially, the sample is completely saturated by the WP. Saturation levels inside the porous sample are controlled by imposing the velocity of pistons that displace the reservoir fluids relative to . As a result of choosing constant, saturation rates remain constant during drainage. Simulations are halted as soon as the NWP penetrates into the WP reservoir, that is, at breakthrough.