ECMOR XI - 11th European Conference on the Mathematics of Oil Recovery

In this paper we present an abstract analysis framework for nonconforming approximations of heterogeneous and anisotropic elliptic operators on general 2D and 3D meshes. The analysis applies to discontinuous Galerkin methods, to the popular MPFA O-method, as well as to hybrid finite volume methods. A number of examples is provided in the paper. The guidelines of the analysis can be summarized as follows: (i) each method is re-written in weak formulation by introducing two discrete gradient operators. Penalty terms may be required for non $H^1$-conformal spaces to ensure coercivity; (ii) a discrete $H^1$-norm is introduced such that a discrete version of the Rellich theorem holds; (iii) an a priori estimate on the discrete solution in the discrete $H^1$-norm is derived using the coercivity of the bilinear form; (iv) the convergence of the discrete problem to the continuous one is eventually deduced, thus proving the convergence of the approximation. In order for the above analysis to hold, the gradient operator applied to the discrete solution should be consistent, whereas the one applied to the test function should be weakly convergent in $L^2$. According to the method, further assumptions on the mesh may be required. This analysis framework easily extends to non-linear problems like the ones encountered in the study of multi-phase Darcy flows, and it allows to weaken regularity assumptions on the coefficients of the problem as well as on the exact solution. From a practical viewpoint, its main interest lies in the possibility to analyze and to design new methods as well as to improve existing ones.

We propose a new single-phase local transmissibility upscaling method for adapted grids in 3-D domains that uses spatially varying and compact multi-point flux approximations (MPFA), based on the VCMP method previously introduced for 2-D Cartesian grids. For each cell face in the coarse upscaled grid, we create a local fine grid region surrounding the face on which we solve three generic local flow problems. The multi-point stencils used to calculate the fluxes across coarse grid cell faces involve up to ten neighboring pressure values. They are required to honor the three generic flow problems as closely as possible while maximizing compactness and ensuring that the flux approximation is as close as possible to being two-point. The resulting MPFA stencil is spatially varying and reduces to two-point expressions in cases without full-tensor anisotropy. Numerical tests show that the method significantly improves upscaling accuracy as compared to commonly used local methods and also compares favorably with a local-global upscaling method. We also present a corrector method that adapts the stencils locally to guarantee that the resulting pressure matrix is an M-matrix. This corrector method is needed primarily for cases where strong permeability contrasts are mis-aligned with the grid locally. The corrector method has negligible impact on the flow accuracy. Finally, we show how the computed MPFA can be used to guide adaptivity of the grid, thus allowing rapid, automatic generation of grids that can account for difficult geologic features such as faults and fractures and efficiently resolve fine-scale features such as narrow channels.

For simulation of three-phase flow processes, such as water-alternating-gas injection (WAG), accurate descriptions of the three-phase capillary pressures and relative permeabilities as functions of the phase saturations are essential. Instead of deriving these functions from experiments, which is very difficult, they may be obtained from pore-scale network simulations. During multiple injection processes, especially under mixed-wet conditions, isolated clusters of the three phases occur. These clusters may still be mobilized, even under the assumption of capillary-dominated flow, as part of so-called multiple displacement chains of adjacent clusters stretching from inlet to outlet of the model. To determine the most favourable displacement chain requires implementation of an efficient shortest path algorithm. For the present problem the distances or costs are the capillary entry pressures between the various phase clusters. Because these entry pressures can also be negative, so-called negative cost cycles arise, which invalidate traditional algorithms. Instead, we have implemented an efficient shortest path algorithm with negative cycle detection If negative cycles arise, the corresponding cyclic displacement chains are carried out before any linear displacement from inlet to outlet. Negative cost cycles correspond to spontaneous displacements, for which the prevailing invading phase pressure is higher than required. Additionally, and probably even more significant, the shortest path algorithm determines in principle the pressures of all phase clusters in the network. These pressures are in turn used in the accurate calculation of the three-phase entry pressures, as well as the film and layer volumes and conductances in pore corners. Some WAG network simulations in mixed-wet media have been carried out using the new algorithm to demonstrate the occurrence of long displacement chains and the much improved efficiency, also for increasing network sizes. A significant number of negative cycles occurred, but the corresponding displacement chains were short. Additionally, we have analysed the obtained values of the cluster pressures.

Fractional-flow theory has proven useful for understanding the factors that control foam displacements and as a benchmark against which to test foam simulators. Most applications of fractional-methods to foam have excluded oil. Recently, Mayberry and Kam (SPE 100964) presented out fractional-flow solutions for foam injection with a constant effect of oil on the foam. We extend fractional-flow methods to foam displacements with oil, using the effects of oil and water saturations on foam as represented in the STARS simulator. There can be abrupt shifts in the composition paths at the limiting water and oil saturations for foam stability. In the immiscible three-phase SAG displacements examined, if foam collapses at the initial oil saturation in the reservoir, there is a very-small-velocity shock from the injected condition to complete foam collapse. The displacement is nearly as inefficient as if no foam were present at all. It does not matter in these cases whether foam is weakened by low water saturation. The displacement is efficient, however, if foam is unaffected by oil but weakened at low water saturation. These results may reflect our foam model, where foam is only partially destroyed at low water saturations but is completely destroyed by high oil saturation. Two idealized models for three-phase displacements can be represented at two-phase displacements with chemical or miscible shocks: A model for a first-contact miscible gas flood with foam suggests an optimal water fraction in foam that puts the gas front just slightly ahead of the foam (surfactant) front. An idealized model of a surfactant flood with foam for mobility control suggests it is important to inject a sufficiently high water fraction in the foam that the gas front is behind the surfactant front as the flood proceeds. We present simulations to verify the solutions obtained with fractional-flow methods.

Growing interest in understanding, predicting, and controlling advanced oil recovery methods emphasizes the importance of numerical methods that exploit the nature of the underlying physics. The Fully Implicit Method offers unconditional stability in the sense of discrete approximations. This stability comes at the expense of transferring the inherent physical stiffness onto the coupled nonlinear residual equations which need to be solved at each time-step. Current reservoir simulators apply safe-guarded variants of Newton’s method, and often can neither guarantee convergence, nor provide estimates of the relation between convergence rate and time-step size. In practice, time-step chops become necessary, and are guided heuristically. With growing complexity, such as in thermally reactive compositional models, this can lead to substantial losses in computational effort, and prohibitively small time-steps. We establish an alternate class of nonlinear iteration that both converges, and associates a time-step to each iteration. Moreover, the linear solution process within each iteration is performed locally. By casting the nonlinear residual for a given time-step as an initial-value-problem, we formulate a solution process that associates a time-step size with each iteration. Subsequently, no iterations are wasted, and a solution is always attainable. Moreover, we show that the rate of progression is as rapid as a standard Newton counterpart whenever it does converge. Finally, by exploiting the local nature of nonlinear waves that is typical to all multiphase problems, we establish a linear solution process that performs computation only where necessary. That is, given a linear convergence tolerance, we identify the minimal subset of solution components that will change by more than the specified tolerance. Using this a priori criterion, each linear step solves a reduced system of equations. Several challenging large-scale simulation examples are presented, and the results demonstrate the robustness of the proposed method as well as its performance.

The analysis of the multi point flux approximation (MPFA) method has so far relied on the possibility of seeing it as a mixed finite element method for which the convergence is then established. This type of analysis has been successfully applied to triangles and quadrilaterals, lately also in the case of rough meshes. Another well known conservative method, the mimetic finite difference method, has also traditionally relied on the analogy with a mixed finite element method to establish convergence. Recently however a new type of analysis proposed by Brezzi, Lipnikov and Shashkov, cf. [3], permits to show convergence of the mimetic method on a general polyhedral mesh. We propose to formulate the MPFA O-method in a mimetic finite difference framework, in order to extend the proof of convergence to polyhedral meshes. The general nonsymmetry of the MPFA method and its local explicit flux, found trough a splitting of the flux, are challenges that require special attention. We discuss various qualities needed for a structured analysis of MPFA, and the effect of the lost symmetry.

Streamline simulation is particularly efficient for incompressible two-phase flow, for which one can use front tracking to solve the 1-D transport problems along streamlines. Here we investigate the extension of this method to compressible (and immiscible) flow and discuss some of the difficulties involved, and in particular the choices one has in writing the 1-D transport equation(s). Our study is motivated by the simulation of CO2 injection, and we therefore also develop methods that are particularly suited for solving compressible flow where one phase is incompressible. Altogether, we present four front-tracking methods that are based on a combination of solving ordinary Riemann problems and Riemann problems with discontinuous flux.

A key aspect in any algebraic multilevel procedure is to be able to reliably capture the physical behavior behind system coefficients. However, the modeling of complex reservoir scenarios generally involves the computation of extremely discontinuous and nonlinear coefficients that, in turn, compromise the robustness and efficiency of smoothing and coarsening strategies. In addition, the need of dealing with large discretized domains leads to highly ill-conditioned systems that represent a true challenge to any linear solver technology known today. In this work, we exploit the fact that flow trend information can be used as a basis for developing a robust percolative aggregation (PA) two-stage preconditioning method. To this end, we identify coefficient aggregates by a means of an efficient version of the Hoshen-Kopelman percolation algorithm suitable for any flow network structure. By partitioning and reordering unknowns according to these aggregates, we obtain a set of high-conductive blocks followed by a low-conductive block. Diagonal scaling allows for weakening the interaction between these high- and low- conductive blocks plus improving the overall conditioning of the system. The solution of the high-conductive blocks approximates well the solution of the original problem. Remaining sources of errors from this approximation are mainly due to small eigenvalues that are properly eliminated with a deflation step. The combination of the algebraic solution of the high-conductive blocks with deflation is realized as a two-stage preconditioner strategy. The deflation stage is carried out by further isolating the aggregate blocks with a matrix compensation procedure. We validate the performance of the PA approach against ILU preconditioning. Preliminary numerical results indicate that the PA two-stage preconditioning can be used as a promising alternative to employ existing algebraic multilevel methods in a more effective way.

Control-volume discretizations for fluid flow in anisotropic porous media are investigated. The method uses numerical solutions of the simplest model equation, –div[K(x) grad p(x)] = f(x). The permeability tensor, K(x), is allowed to have discontinuities. Multipoint Flux Approximations (MPFA) are used, and transmissibility coefficients are obtained, in the usual way, from local numerical flow experiments (transmissibility upscaling) for each cell face. For regular K-orthogonal grids, with a uniform permeability tensor, the scheme reduces to the standard two-point flux approximation. Monotonicity of the solution matrix is discussed and a version of the method that provides an M-matrix is described. This discretization scheme is applied to reservoir simulation on 3D structured grids with distorted geometry, highly anisotropic media, and discontinuities in the permeability tensor. Simulation results are presented and compared with results from other two-point and multipoint flux approximations.

Upscaling is commonly applied to account for the effects of fine-scale permeability heterogeneity in coarse-scale simulation models. Ensemble level upscaling (EnLU) is a recently developed approach that aims to efficiently generate coarse-scale flow models capable of reproducing the ensemble statistics (e.g., P50, P10 and P90) of fine-scale flow predictions for multiple reservoir models. Often the most expensive part of standard coarsening procedures is the generation of upscaled two-phase flow functions (e.g., relative permeabilities). EnLU provides a means for efficiently generating these upscaled functions using stochastic simulation. This involves the use of coarse-block attributes that are both fast to compute and representative of the effects of fine-scale permeability on upscaled functions. In this paper, we establish improved attributes for use in EnLU, namely the coefficient of variation of the fine-scale single-phase velocity field (computed during computation of upscaled absolute permeability) and the integral range of the fine-scale permeability variogram. Geostatistical simulation methods, which account for spatial correlations of the statistically generated upscaled functions, are also applied. The overall methodology thus enables the efficient generation of coarse-scale flow models. The procedure is tested on 3D well-driven flow problems with different (Gaussian) permeability distributions and high fluid mobility ratios. EnLU is shown to capture the ensemble statistics of fine-scale flow results (flow rate and oil cut as a function of time) with similar accuracy to full flow-based methods but at a small fraction of the computational cost.