An Efficient Solver for Nonlinear Multiphase Flow Based on Adaptive Coupling of Flow and Transport

Fully Implicit discretization of flow and transport equations gives rise to a system of coupled nonlinear equations that is typically solved using standard Newton method. For a given timestep size, even if the Newton-based iterative procedure converges, the cost associated with updating all the unknowns simultaneously can be quite expensive. Conventional sequential-implicit strategies can be used to reduce the cost, but they suffer from severe restrictions on the allowable timestep size.

In this paper we formulate, verify and analyze the computational efficiency of a new nonlinear solution strategy. The crux of the proposed algorithm is the use of a hybrid strategy to treat the co-current and counter-current flow regions differently. At each Newton iteration, we first update the pressure variables by solving the Schur complement reduced form of the equations. To avoid costly computation of the reduced Jacobian, we employ phase-based potential ordering, giving rise to a lower triangular Jacobian that can be used to compute the reduced Jacobian efficiently. After updating pressures, we decompose the domain into components in a way that can be sorted and traversed from upstream to downstream. A component in the co-current flow regions is made up of a single cell, whereas a component made up of multiple connected cells indicates counter-current flow in which pressure and saturation are tightly coupled. Marching down from upstream to downstream, for a single-cell component we update saturation nonlinearly by solving the scalar saturation equation(s). For a multi-cell component we discard the pressure of the corresponding cells obtained in the first step, and then perform a simultaneous linear update of the saturation and pressure variables by solving the local linear system. Once all of the components are visited and updated, the iteration is over.

We present a variety of challenging numerical examples in 2-D and 3-D in the presence of strong gravity and heterogeneity. Our results show that as compared with standard Newton method, the proposed hybrid solver has a lower computational cost per iteration without compromising the allowable timestep size. The computational efficiency gain depends on the number and size of the components which vary over the course of iterations. At best, pressure and saturation updates are fully decoupled and saturation variables are updated sequentially one at a time. We present a rigorous complexity analysis of the algorithm for linear solvers with arbitrary order of complexity.