1. Introduction
Fluid-fluid interfacial area plays an important role in many practical applications related to two-phase flow and mass transfer, such as enhanced oil recovery, CO2 sequestration, and groundwater contamination by non-aqueous phase liquids (NAPLs). This repository presents a reference workflow to investigate the influence of displacement regime, geometric disorder of the porous medium, and wettability on the fluid-fluid interfacial area between two phases.
2 Geometry generation
The computational domain is a two-dimensional rectangular porous medium (with a length of 60 mm and a width of 20 mm), constructed by placing 1000 non-overlapping spheres as the solid phase on a triangular lattice with an identical lattice spacing of d = 1.386 mm. To investigate the effect of geometric disorder on the fluid-fluid interfacial area between two phases, λ is introduced by assigning the radius of each sphere randomly from a uniform distribution, r ∈ [1 - λ, 1 + λ]r̄, where r̄ is the mean sphere radius.
Using this method, four levels of geometric disorder, namely λ = 0, 0.1, 0.2, and 0.3, are generated. For all porous media considered, the porosity is maintained constant at φ = 0.5 to isolate the effect of geometric disorder. These geometries can be found in the folder porous_media_geometries. To characterize the pore space influenced by the level of geometric disorder, a pore network model (PNM) is extracted from the porous medium geometry using PoreSpy. Specifically, the SNOW (Subnetwork of the Oversegmented Watershed) algorithm is employed to identify individual pores and throats, providing key parameters such as pore and throat radii and connectivity.
3. Numerical implementation
The typical Volume of Fluid (VOF) hydrodynamics solver, interFoam, within the OpenFOAM (ESI, v2306) framework, is employed to simulate two-phase fluid flow in porous media. The pressure-velocity coupling is handled using the extended PIMPLE algorithm, which combines the iterative pressure-correction approach of SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) with the transient acceleration of PISO (Pressure-Implicit with Splitting of Operators).
It should be noted that an adaptive time-stepping strategy is recommended to ensure both numerical stability and computational efficiency. In this approach, the time step is dynamically adjusted during the simulation to maintain the Courant number below 0.2 throughout the computational domain. An immiscible two-phase flow system is considered, where water serves as the defending wetting phase (WP) and the other fluid as the invading non-wetting phase (NWP). The solid surface is assumed to be water-wet. The interfacial tension between the two fluids is σ = 0.05 N/m, and the density of the NWP is 703 kg/m3.
In order to study the effect of the displacement regime, the viscosity of the invading NWP is varied to achieve different combinations of Capillary number and the viscosity ratio between the invading NWP and the defending WP. The viscosity of the invading NWP can be modified in constant/transportProperties in the provided OpenFOAM example case (OpenFOAM_example_case). Similarly, the effect of wettability is investigated by prescribing different contact angles in 0/alpha.water.
4. Post-processing
Using the VOF method, the invading NWP and WP are identified based on the volume fraction α, with the interfacial area defined by the iso-contour α = 0.5. The interfacial area is quantified by extracting the length of the α = 0.5 iso-contour in the two-dimensional domain.
All post-processing operations, including the connectivity analysis of WP clusters and the extraction of interfacial areas, are performed using the open-source Python package PyVista. The post-processing codes can be found in the folder postProcessing. Before running the post-processing scripts, the OpenFOAM simulation results should be converted into VTK format using the foamToVTK command within the OpenFOAM environment.
