Influence of Fractal Dimension on Gas-Driven Two-Phase Flow in Fractal Porous Media: A VOF Model-Based Simulation

Gas-liquid two-phase flow in fractal porous media is pivotal for engineering applications, yet it remains challenging to be accurately characterized due to complex microstructure-flow interactions. This study establishes a pore-scale numerical framework integrating Monte Carlo-generated fractal porous media with Volume of Fluid (VOF) simulations to unravel the coupling among pore distribution characterized by fractal dimension (Df), flow dynamics, and displacement efficiency. A pore-scale model based on the computed tomography (CT) microstructure of Berea sandstone is established, and the simulation results are compared with experimental data. Good agreement is found in phase distribution, breakthrough behavior, and flow path morphology, confirming the reliability of the numerical simulation method. Ten fractal porous media models with Df ranging from 1.25~1.7 were constructed using a Monte-Carlo approach. The gas-liquid two-phase flow dynamics was characterized using the VOF solver across gas injection rates of 0.05–5 m/s, in which the time-resolved two-phase distribution patterns were systematically recorded. The results reveal that smaller fractal dimensions (Df = 1.25~1.45) accelerate fingering breakthrough (peak velocity is 1.73 m/s at Df = 1.45) due to a bimodal pore size distribution dominated by narrow channels. Increasing Df amplifies vorticity generation by about 3 times (eddy viscosity is 0.033 Pa·s at Df = 1.7) through reduced interfacial curvature, while tortuosity-driven pressure differentials transition from sharp increases (0.4~6.3 Pa at Df = 1.25~1.3) to inertial plateaus (4.8 Pa at Df = 1.7). A nonlinear increase in equilibrium gas volume fraction (fav = 0.692 at Df = 1.7) emerges from residual gas saturation and turbulence-enhanced dispersion. This behavior is further modulated by flow velocity, with fav peaking at 0.72 under capillary-dominated conditions (0.05 m/s), but decreasing to 0.65 in the inertial regime (0.5 m/s). The work quantitatively links fractal topology to multiphase flow regimes, demonstrating the critical role of Df in governing preferential pathways, energy dissipation, and phase distribution.