An Efficient Fully-Implicit MFD-MUSCL Method Based on a Novel Multislope Limiting Procedure

Abstract Standard reservoir simulation schemes use single-point upstream weighting for computing the convective fluxes when multi-phase/component fluids are present. These schemes are only first-order and may cause high viscosity effect. Second-order schemes provide a better resolution and reduce the smoothing near the shocks. In reservoir simulation practice, implicit discretisations capable of taking large time steps are preferred over explicit schemes. However, even for a simple first-order method, solving a large nonlinear system is often very expensive; Extra coupling and nonlinearity of the discretised equations can be introduced from a higher-order spatial discretisation. It has been shown that strong nonlinearity as well as lack of continuous differentiability in numerical flux function and flux limiter can cause serious nonlinear convergence problem. Cell-centered finite-volume (CCFV) discretizations may offer several attractive features, especially for fluid flow in heterogeneous or fractured porous media. The objectives of this work are to develop a fully-implicit CCFV framework that could achieve second-order spatial accuracy on smooth solutions, while at the same time maintain robustness and nonlinear convergence performance. We develop a novel multislope MUSCL method which interpolates the required values at the edge centroids in a simpler way by taking advantage of some geometric features of the triangular mesh. Through the edge centroids, the numerical diffusion caused by mesh skewness is significantly reduced, and optimal second-order accuracy can be achieved. An improved gradually-switching piecewise-linear flux-limiter is introduced according to mesh non-uniformity in order to prevent spurious oscillations. The smooth flux-limiter can achieve high accuracy without degrading nonlinear convergence behavior. For the discretization of pressure and Darcy velocities, a mimetic finite difference method that provides flux- continuity and an accurate total velocity field is used. The fully-coupled MFD-MUSCL framework is adapted to accommodate a lower-dimensional discrete fracture-matrix (DFM) model. Several numerical tests with discrete fractured system are carried out to demonstrate the efficiency and robustness of the numerical model. The results show that the high-order MUSCL method effectively reduces numerical diffusion, leading to improved resolution of saturation fronts compared with the first-order method. In addition, it is shown that the developed multislope method and adaptive flux limiter exhibit superior nonlinear convergence compared with other alternatives.