Convergent sequential fully implicit method for reservoir simulation

We propose a stable time stepping algorithm for sequential fully implicit method (SFIM) in solving the coupled flow and transport equations for multi-phase flow in reservoir simulation. In developing a stable sequential fully implicit method, we are concerned with the flow and transport systems that yield a stable, unique solution for a sufficiently small time step size. The sequential fully implicit method comprises three major steps: (1) pressure solution, (2) conservative total velocity calculation, and (3) saturation solution from the nonlinear transport equation. These three major steps are coupled by outer-loop iterations. In solving the transport equation, the fractional flow formulation with fixed total velocity is employed. We assume that the transport equation can be always solved with an underrelaxation method for a given conservative velocity field. For example, by analyzing the S-shaped fractional flow curve, Jenny et al. (2009) [17] developed an unconditionally stable under-relaxation method for the iterative solution of two-phase flow with weak gravity effect. In multi-phase (two or three phases) flow with gravity effect, the fluxes between two cells are split into viscous and gravitational fluxes. The fractional flow for viscous flux is a function of the upwind cell saturation, whereas that for gravitational flux is a function of two saturations of the neighboring cells. For a simple representation of multivariable functions in Newton iterations, the fractional flow curves are parameterized along the Newton iterative direction. Since the fractional flow curves for viscous and gravitational flows are distinctly different, we devise two under-relaxation methods for multi-phase flow with gravity effect. In SFIM, the algorithm can be optimized based on the total velocity changes during a time step. In general, one iteration of saturation is computationally less expensive than that for one pressure calculation. As a result, if the total velocity is not a strong function of saturations (e.g., viscous flow with small gravity number), for one outer-loop calculation, the transport equation is iterated to convergence. When the total velocity is a strong function of saturations (e.g., lock-exchange problem), saturations are computed only for one iterative update for each outer-loop iteration. It ensures that the solution path will closely follow the trust region of pressure and saturation in Newton's method. We demonstrate convergence of SFIM with numerical examples: gravity segregation, a linear displacement model, reservoir models with wells, a lock-exchange problem, etc. (c) 2023 Elsevier Inc. All rights reserved.