Multiscale Finite-Volume CVD-MPFA Formulations on Structured and Unstructured Grids

This paper presents the development of finite-volume multiscale methods for quadrilateral and triangular unstructured grids. Families of Darcy-flux approximations have been developed for consistent approximation of the general tensor pressure equation arising from Darcy's law together with mass conservation. The schemes are control-volume distributed (CVD) with flow variables and rock properties sharing the same control-volume location and are comprised of a multipoint flux family formulation (CVD-MPFA). The schemes are used to develop a CVD-MPFA based multiscale finite-volume (MSFV) formulation applicable to both structured and unstructured grids in two dimensions. The basis functions are a key component of the MSFV method, and are a set of local solutions, usually defined subject to Dirichlet boundary conditions. A generalization of the Cartesian grid Dirichlet basis functions described in [P. Jenny, S. H. Lee, and H. A. Tchelepi, J. Comput. Phys., 187 (2003), pp. 47--67] is presented here for unstructured grids. Whilst the transition from a Cartesian grid to an unstructured grid is largely successful, use of Dirichlet basis functions can still lead to pressure fields that exhibit spurious oscillations in areas of strong heterogeneity. New basis functions are proposed in an attempt to improve the pressure field solutions where Neumann boundary conditions are imposed almost everywhere, except corners which remain specified by Dirichlet values

Interior boundary-aligned unstructured grid generation and cell-centered versus vertex-centered CVD-MPFA performance

Grid generation for reservoir simulation must honor classical key constraints and ensure boundary alignment such that control-volume boundaries are aligned with geological features including layers, shale barriers, fractures, faults, pinch-outs, and multilateral wells. Novel unstructured grid generation methods are proposed that automate control-volume and/or control point boundary alignment and yield perpendicular-bisector (PEBI) meshes both with respect to primal and dual (essentially PEBI) cells. In order to honor geological features in the primal configuration, we introduce the idea of protection circles that contain segments of key geological boundaries, while in order to generate a dual-cell feature aligned grid, we construct halos around key geological features. The grids generated are employed to study comparative performance of cell-centred versus cell-vertex flux-continuous control-volume distributed multi-point flux approximation (CVD-MPFA) finite-volume formulations using equivalent degrees of freedom and thus ensure application of the most efficient methods. The CVD-MPFA formulation (c.f. Edwards et al.) in cell-centred and cell-vertex modes is somewhat analogous and requires switching control-volume from primal to dual or vice versa, together with appropriate data structures and boundary conditions, however dual-cells are generated after primal grid generation. The relative benefits of both types of approximation, i.e., cell-centred versus vertex-centred, are contrasted in terms of flow resolution and degrees of freedom required

Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model

A novel cell-centered control-volume distributed multi-point flux approximation (CVD-MPFA) finite-volume formulation is presented for discrete fracture-(rock)matrix flow simulations. The grid is aligned with the fractures and barriers which are then modeled by lower-dimensional interfaces located between rock matrix cells in the physical domain. The n D (n-dimension) pressure equation in the rock matrix is coupled with the (n−1)D pressure equation in the fractures, leading to the first reduced dimensional flux-continuous CVD-MPFA formulation. This formulation naturally handles fractures efficiently on unstructured grids. Matrix-fracture fluxes are expressed in terms of matrix and fracture pressures, resulting in a transfer function, which is added to the lower-dimensional flow equation. An additional transmission condition is used between matrix cells separated by low permeable fractures to couple the velocity and pressure jump across the fractures. Numerical tests serve to assess the convergence and accuracy of the lower-dimensional fracture model for lower anisotropic fractures having a range of apertures and permeability tensors. A tracer flow transport equation is solved for problems with single and intersecting fractures. A lower-dimensional mass balance for intersecting fracture cells circumvents the more restrictive CFL condition corresponding to standard equi-dimensional approximation with explicit time discretization. Lower-dimensional fracture model results are compared with hybrid-grid and equi-dimensional model results. Fractures and barriers are efficiently modeled by lower-dimensional interfaces which yield comparable results to those of the equi-dimensional model. Highly conductive fractures are modeled as lower-dimensional entities without the use of locally refined grids that are required by the equi-dimensional model, while pressure continuity across fractures is built into the model, without depending on the extra degrees of freedom which must be added locally by the hybrid-grid method. The lower-dimensional fracture model also yields improved results when compared to those of the hybrid-grid model for fractures with low-permeability in the normal direction to the fracture where pressure is discontinuous. In addition, transient pressure simulation involving geologically representative complex fracture networks is presented

Simulation of multi-component multi-phase fluid flow in two-dimensional anisotropic heterogeneous porous media using high-order control volume distributed methods

Peer reviewedPostprin

Quasi-positive families of flux continuous finite volumes schemes in two and three dimensions.

In this thesis, new families of full pressure support flux-continuous, locally conservative, finite-volume schemes are presented for solving the general geometry-permeability tensor pressure equation on structured and unstructured grids in two and three dimensions. The families of flux-continuous schemes have also been referred to in the literature as Multi-point Flux Approximation or MPFA schemes. The schemes are applicable to the general tensor pressure equation with discontinuous coefficients and remove the 0(1) errors introduced by standard reservoir simulation (two-point flux) schemes when applied to full, anisotropic and asymmetric permeability tensor flow approximation. Such tensors may arise when fine scale permeability distributions are upscaled to obtain gridblock-scale permeability distributions. In contrast to the previous MPFA schemes which assume point-wise pressure and flux continuity locally, the new families of schemes presented in the work recover full pressure continuity across the interface between neighboring subcells. The M-matrix conditions [1, 2] define the upper limits for ensuring a local maximum principle is obtained for full-tensor fields. A key condition is that the modulus of the off-diagonal tensor coefficients are bounded by the minimum of the diagonal coefficients. For higher anisotropic ratios, when the resulting discrete matrices violate these bounds these schemes can violate the maximum principle (as with more standard methods) and the numerical pressure solutions can consequently exhibit spurious oscillations. The new family of schemes yield improved performance for challenging problems where earlier flux-continuous schemes exhibit strong spurious oscillations. The M- matrix analysis leads to an optimal quadrature range for these methods. The degree of freedom within the family of full pressure continuity schemes presented is shown to maximise the quadrature range of the flux-continuous schemes. For strongly anisotropic full-tensor cases where M-matrix conditions are violated, it is shown that the earlier families of schemes cannot avoid decoupling of the solution which leads to severe spurious oscillations in the discrete solution. The full quadrature range of the new schemes permits use of quadrature points that were previously out of range for the earlier methods, and that the resulting schemes minimize spurious oscillations in discrete pressure solutions. The new formulation leads to a more robust quasi-positive family of flux-continuous schemes applicable to general discontinuous full-tensor fields. This work also extends the single parameter family of FPS schemes to double families of schemes with general flexibility in quadrature that allow different quadrature points to be used on different control-volume subfaces. The new schemes minimize spurious oscillations in discrete pressure solutions. The new formulation leads to more robust quasi-positive families of flux-continuous schemes applicable to general discontinuous full-tensor fields. The full pressure support flux continuous schemes also extend to 3D on structured and unstructured grids. Surface auxiliary control volume and volume auxiliary control volume are introduced to handle extra degrees of freedom which are required for full pressure continuity over neighboring subcell surface. The new schemes are shown to be beneficial in high anisotropic test cases while remaining comparable with previous tetrahedral pressure support (TPS) schemes in terms of convergence rate. Multi-family schemes in 3D are also presented in this work. This is the extension of 2D double family to 3D. Compared to single family FPS schemes, multi-family schemes are shown to be able to maximize the quadrature and have incomparable flexibility over previous schemes, leading to improved solutions

Three-dimensional unstructured gridding for complex wells and geological features in subsurface reservoirs, with CVD-MPFA discretization performance

Grid generation for reservoir simulation, must honour classical key geological features and multilateral wells. The features to be honoured are classified into two groups; (1) involving layers, faults, pinchouts and fractures, and (2) involving well distributions. In the former, control-volume boundary aligned grids (BAGs) are required, while in the latter, control-point (defined as the centroid of the control-volume) well aligned grids (WAGs) are required. Depending on discretization method type and formulation, a choice of control-point and control-volume type is made, i.e. for a cell-centered method the primal grid cells act as control-volumes, otherwise for a vertex-centered method the dual-grid cells act as control-volumes. Novel three-dimensional unstructured grid generation methods are proposed that automate control-volume boundary alignment to geological features and control point alignment to complex wells, yielding essentially perpendicular bisector (PEBI) meshes either with respect to primal or dual-cells depending on grid type. Both grid types use tetrahedra, pyramids, prisms and hexahedra as grid elements. Primal-cell feature aligned grids are generated using special boundary surface protection techniques together with constrained cell-centered well trajectory alignment. Dual-cell feature aligned grids are generated from underlying primal-meshes, whereby features are protected such that dual-cell control-volume faces are aligned with interior feature boundaries, together with protected vertex-centered (control point) well trajectory alignment. The novel methods of grid generation presented enable practical application of both method types in 3-D for the first time. The primal and dual grids generated here demonstrate the gridding methods, and enable the first comparative performance study of cell-vertex versus cell-centered control-volume distributed multi-point flux approximation (CVD-MPFA) finite-volume formulations using equivalent mesh resolution on challenging problems in 3-D. Pressure fields computed by the cell-centered and vertex-centered CVD-MPFA schemes are compared and contrasted relative to the respective degrees of freedom employed, and demonstrate the relative benefits of each approximation type. Stability limits of the methods are also explored. For a given mesh the cell-vertex method uses approximately a fifth of the unknowns used by a cell-centered method and proves to be the most beneficial with respect to accuracy and efficiency. Numerical results show that vertex-centered CVD-MPFA methods outperform cell-centered CVD-MPFA method

CVD-MPFA full pressure support, coupled unstructured discrete fracture–matrix Darcy-flux approximations

Two novel control-volume methods are presented for flow in fractured porous media, involving coupling the control-volume distributed multi-point flux approximation (CVD-MPFA (c.f. Edwards et al.)) constructed with full pressure support (FPS), to two types of discrete fracture-matrix approximation for flow simulation on unstructured grids; (i) involving hybrid grids and (ii) a lower dimensional fracture model. Flow is governed by Darcy's law together with mass conservation both in the rock matrix and in fractures, where large discontinuous permeability tensors can occur. Finite-volume FPS schemes are more robust than the earlier CVD-MPFA triangular pressure support (TPS) schemes for problems involving strongly anisotropic homogeneous and heterogeneous full-tensor permeability fields. We use a cell-centred hybrid-grid method, where fractures are represented by lower-dimensional interfaces between matrix grid cells in the physical mesh, and expanded to equi-dimensional cells in the computational domain. We present a simple procedure to form a consistent hybrid-grid locally for a dual-cell. We also propose a novel hybrid-grid for intersecting fractures, for the FPS method, which improves the condition number of the global linear system and permits larger time steps for tracer transport. The tracer flow transport equation is coupled with the pressure equation and the results provide flow parameter assessment of the fracture models. Transport results obtained via TPS and FPS hybrid-grid formulations are compared with corresponding results of fine-scale explicit equi-dimensional formulations. The results show that the hybrid-grid FPS method applies to general full-tensor fields and provides improved robust approximations compared to the hybrid-grid TPS method for fractured domains, for both weakly anisotropic permeability fields and in particular for very strong anisotropic full-tensor permeability fields where the TPS scheme exhibits spurious oscillations. The hybrid-grid FPS formulation is extended to compressible flow and the results demonstrate the method is also robust for transient flow. Furthermore, FPS is coupled with a lower-dimensional fracture model, where fractures are strictly lower-dimensional in the physical mesh. Comparisons of the hybrid-grid FPS method and the FPS lower-dimensional fracture model are presented for several cases of isotropic and strongly anisotropic fractured media which illustrate the benefits of the respective methods

Multigrid solvers for multipoint flux approximations of the Darcy problem on rough quadrilateral grids

In this work, an efficient blackbox-type multigrid method is proposed for solving multipoint flux approximations of the Darcy problem on logically rectangular grids. The approach is based on a cell-centered multigrid algorithm, which combines a piecewise constant interpolation and the restriction operator by Wesseling/Khalil with a line-wise relaxation procedure. A local Fourier analysis is performed for the case of a Cartesian uniform grid. The method shows a robust convergence for different full tensor coefficient problems and several rough quadrilateral grids.Francisco J. Gaspar has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska–Curie grant agreement no. 705402, POROSOS. The work of Laura Portero is supported by the Spanish project MTM2016-75139-R (AEI/FEDER, UE) and the Young Researchers Programme 2018 from the Public University of Navarre. Andrés Arrarás acknowledges support from the Spanish project PGC2018-099536-A-I00 (MCIU/AEI/FEDER, UE) and the Young Researchers Programme 2018 from the Public University of Navarre. The work of Carmen Rodrigo is supported by the Spanish project PGC2018-099536-A-I00 (MCIU/AEI/FEDER, UE) and the DGA (Grupo de referencia APEDIF, ref. E24_17R)

A Finite-Volume-Based Module for Unsaturated Poroelasticity

In this chapter, we present fv-unsat, a multipoint finite-volume–based solver for unsaturated flow in deformable and nondeformable porous media. The latter is described using the mixed form of Richards’ equation, whereas the former by the equations of unsaturated poroelasticity. The module aims at flexibility, relying heavily on discrete operators and equations, exploiting the automatic differentiation framework provided by the MATLAB Reservoir Simulation Toolbox (MRST). Our examples cover two numerical convergence tests and two three-dimensional practical applications, including the water infiltration process in a nondeformable soil column and a realistic desiccation process of a deformable clay sample using atmospheric boundary conditions. The resulting convergence rates are in agreement with previously reported rates for single-phase models, and the practical applications capture the physical processes accurately.publishedVersio