A new Timestep Criterion for the Simulation of Immiscible Two-Phase Flow with IMPES Solvers

TL;DR

This work addresses stability and efficiency in iterative IMPES simulations of immiscible two-phase flow with compressible phases by introducing a generalized characteristic wave velocity timestep criterion that leverages numerically computed derivatives of the wetting-velocity with respect to saturation. The FV-based solver combines a total-velocity discretization, consistent saturation- and pressure-updates, and a CFL-bound derived from $\omega^{k,max}_{ij}$ to adapt $\Delta t$ across capillary-dominated, gravity-driven, and pressure-drop-driven regimes. Across capillary rise, Buckley-Leverett, capillary gravity equalization, gas compression, and discontinuous material cases, the generalized criterion reduces time iterations compared with the Coats criterion while maintaining comparable accuracy, and it improves non-wetting phase mass conservation for compressible flows. The results demonstrate robustness and efficiency without parameter tuning, making the approach broadly applicable to complex porous-media simulations with discontinuities and compressibility.

Abstract

We present an iterative IMPES solver and a novel timestep criterion for the simulation of immiscible two-phase flow involving compressible fluid phases. The novel timestep criterion uses the Courant-Friedrichs-Lewy (CFL) condition and employs numerically computed velocity derivatives to adapt the timestep size, regardless of the dominant flow characteristics. The solver combined with this timestep criterion demonstrates both efficiency and robustness across a range of flow scenarios, including pressure drop dominated and capillary dominated flows with compressible and incompressible fluid phases, without the need to adjust any numerical parameters. Furthermore, it successfully reaches the expected stationary states in a case involving discontinuous porous media parameters such as porosity, permeabilities, and capillary pressure function. Comparison with the established Coats timestep criterion reveals that our approach requires fewer time iterations while maintaining comparable accuracy on the Buckley-Leverett problem and a gravity-capillary equalization example with a known stationary state. Additionally, in an example with air compression, the new timestep criterion leads to a significantly improved non-wetting phase mass conservation compared to the Coats criterion.

Figures (14)

• Figure 1: This is a cutout of a two-dimensional grid with the relevant grid parameters as defined in Section \ref{['section_numerics']}.
• Figure 2: The image displays the saturation profiles from the capillary rise simulations at $10000 s$. The results were obtained using the characteristic wave velocity criterion, the generalized characteristic wave velocity criterion, and the Coats criterion. The grid used measures $10 \times 100$ voxels.
• Figure 3: The numerical and analytic solutions of the Buckley-Leverett problem at various times are shown. Solid lines correspond to the analytic solution, while dashed lines indicate the numerical results. The geometry is discretized with a voxel size of $10^{-3} m$, and the simulations are performed using $C_{\text{stab}} = 1$.
• Figure 4: This image shows both the analytic and numerical solutions of the Buckley-Leverett problem at time $t=450 s$. The geometry is discretized with a voxel size of $10^{-3} m$ and the numerical solutions are computed using various stability constants $C_{\text{stab}}$.
• Figure 5: These two images display the errors of the numerical solutions for the Buckley-Leverett problem across different voxel sizes. The stability constant is set to $C_{\text{stab}} = 1$. Black dots indicate the numerical errors, with the left image presenting the $L_1$ error and the right image showing the $L_2$ error. The red lines correspond to linear least-square regression lines, estimating the order of convergence of the solver. These regression lines are computed by a least-squares fit of the ten-logarithm of the voxel sizes and the ten-logarithm of the corresponding error values. The slope of this fit represents the approximate order of convergence, yielding about $0.92$ for the $L_1$ error and $0.48$ for the $L_2$ error.