An Efficient Fully-implicit High-resolution MFD-MUSCL Method for Two-phase Flow with Gravity in Discrete Fractured Media

Standard reservoir simulation schemes employ single-point upstream weighting for approximation of the convective fluxes when multiple phases or components are present. These schemes are only first-order and give a poor approximation and induce high viscosity effect. A second-order scheme provides a better approximation and manages to reduce the viscous smoothing effect in the vicinity of the shocks. In reservoir simulation practice, implicit discretisations capable of taking large time steps are preferred in practical computations. However, assembling and solving a large nonlinear system is often very expensive, even for a simple first-order method, and using a higher-order spatial discretisation introduces extra couplings and increases the nonlinearity of the discretised equations. It has been shown that the strong nonlinearity as well as the 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 the fluid flow in discrete fractured media. The objectives of this work are to develop a novel cell-centered multislope MUSCL method and an adaptive limiting strategy that have improved computational efficiency, smoothness properties, and accuracy. The reconstruction scheme interpolates the required values at the edge centroids in a more straightforward and effective way by making better use of the triangular mesh geometry. Because optimal second-order accuracy can be reached at the edge centroids, the numerical diffusion caused by mesh skewness is also significantly reduced. An improved gradually-switching piecewise-linear flux-limiter is introduced according to mesh non-uniformity in order to prevent spurious oscillations. The developed 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 developed fully-coupled MFD-MUSCL CCFV framework is adapted to accommodate a lower-dimensional discrete fracture-matrix 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 scheme and adaptive flux limiter exhibit superior nonlinear convergence compared with other alternatives.