Nonlinear Analysis of Newton-based Solvers for Multiphase Transport in Porous Media with Viscous and Buoyancy Forces

Numerical simulation of multiphase flow in porous geological formations is widely used for subsurface flow management, including oil and gas production. Large-scale simulation is often limited by the computational speed, where nonlinear convergence is one of the main bottlenecks. The nonlinearity of flow properties, the coupling of driving forces for fluid migration, and the heterogeneities of the formation are three main causes for convergence difficulties. Here, we analyze the nonlinearity of two-phase transport in porous media, and we propose an efficient nonlinear solver based on this understanding.

We focus on immiscible, incompressible, two-phase transport in the presence of viscous and buoyancy forces. We investigate the nonlinearity of the discrete transport equation obtained from finite-volume discretization with Single-Point Upstream weighting (SPU), which is the industry standard. In particular, we study the discretized numerical flux expressed as a function of the upstream and downstream saturations of the total velocity. We analyze the locations and complexity of the unit-flux, zero-flux, and inflection lines of the numerical-flux saturation space. The unit- and zero-flux lines, referred to as kinks, correspond to a switch in the flow directions of the different phases, and if SPU is used, then the numerical flux is not differentiable at those points. These kinks and inflection lines are major sources of nonlinear convergence difficulties, especially when their locations in the numerical flux depend on both saturations in the upstream and downstream cells. Our analysis of the discrete (numerical) flux offers a theoretical basis of the convergence challenges associated with multi-cell problems and serves as a foundation for developing efficient nonlinear solvers.

With this understanding, we propose a nonlinear solution scheme that is a significant refinement of the works of Jenny et al. (2009) and Wang and Tchelepi (2013) . We divide the flux function into saturation ‘trust regions’ delineated by the kinks and inflection lines. Determining the boundary of each trust region is straightforward, and it only needs to be computed once in a preprocessing step before performing a simulation. The Newton updates are performed such that two successive iterations do not cross any trust-region boundary. If a crossing is detected, the saturation value is chopped back to the boundary. Our saturation chopping algorithms captures the inflection lines of the numerical flux much more accurately than the treatment of Wang and Tchelepi (2013) . The nonlinear convergence behavior is analyzed using numerical examples, and significant improvements over existing trust-region nonlinear solvers are demonstrated.