Primal Grid Research Articles

A very high-order flux reconstruction approach coupled to the MPFA-QL finite volume method for the numerical simulation of oil-water flows in 2D petroleum reservoirs

In recent years, high-order methods for the discretization of the advection dominated saturation equation, and that results from the modelling of the 2D Oil-Water displacement through porous formations, have received renewed attention from the CFD community, as a suitable alternative to the traditional first-order upwind (FOU) finite volume (FV) scheme, which is the standard in industry. A key aspect of high-order methods applied to numerical simulation of petroleum reservoirs is the lower computational cost for a comparable degree of accuracy against their lower-order counterparts. In this context, in the present paper, for the first time in literature, we propose the use of the very high order Flux Reconstruction (FR) formulation for the discretization of the advective terms in the saturation equation using unstructured quadrilateral meshes. To solve the elliptic pressure equation, we use the recently developed Multipoint Flux Approximation method with a Quasi-Local Stencil (MPFA-QL). This method is very robust and able to handle highly heterogeneous and anisotropic permeability tensors. The system of equations is solved in a segregated manner by using the classical Implicit Pressure Explicit Saturation (IMPES) procedure. To properly couple the cell-centered MPFA-QL method with the FR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. As the cell-centered finite volume method naturally delivers fluxes across control surfaces of the primal grid, then, we use a reconstruction operator based on the first-order Brezzi–Douglas–Marini (BDM) space, and the Piola transformation, to get the complete knowledge of the velocity field throughout the domain. To eliminate spurious oscillations around the discontinuities, which are typical in high-order methods, we have used a hierarchical multidimensional limiting process (hMLP). To compute fluxes over the control surfaces we have used the Roe-E approximate Riemann solver. Finally, through some representative examples, we show the accuracy, efficiency and robustness of our new high order numerical methodology to model oil-water displacements in highly heterogeneous and anisotropic petroleum reservoirs.

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.

Numerical simulation of two-phase flows in 2-D petroleum reservoirs using a very high-order CPR method coupled to the MPFA-D finite volume scheme

The focus of the present work is to model the two-phase flow of oil and water in 2-D anisotropic petroleum reservoirs using a novel appropriate combination of a non-orthodox cell-centered Multipoint Flux Approximation with a Diamond Stencil (MPFA-D) finite volume method to discretize the pressure equation coupled to the very high-resolution Correction Procedure via Reconstruction (CPR) scheme for the discretization of the saturation equation. To suppress numerical oscillations (under/overshoots) near shocks, that are typical in higher-order schemes, and to deliver high accuracy in smooth regions of the solution, a hierarchical Multidimensional Limiting Process (MLP) is used in the reconstruction stage. The integration in time is carried out using a third-order Runge–Kutta method. The coupling of the pressure-saturation set of equations is carried out using a classical Implicit Pressure Explicit Saturation (IMPES) procedure. To properly couple the MPFA-D method with the CPR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. Because the cell-centered finite volume method naturally delivers fluxes across control surfaces of the primal grid, then, we use a reconstruction operator based on the lowest order Raviart–Thomas interpolation functions and the Piola transformation, to get the complete knowledge of the velocity field throughout the domain. Finally, some representative problems were solved to evaluate the accuracy, efficiency, and shock-capturing capabilities of our new numerical methodology.

A pressure-based semi-implicit space–time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier–Stokes equations at all Mach numbers

We propose a new arbitrary high order accurate semi-implicit space–time discontinuous Galerkin (DG) method for the solution of the two and three dimensional compressible Euler and Navier–Stokes equations on staggered unstructured curved meshes. The method is pressure-based and semi-implicit and is able to deal with all Mach number flows.The new DG scheme extends the seminal ideas outlined in [1], where a second order semi-implicit finite volume method for the solution of the compressible Navier–Stokes equations with a general equation of state was introduced on staggered Cartesian grids. Regarding the high order extension we follow [2], where a staggered space–time DG scheme for the incompressible Navier–Stokes equations was presented. In our scheme, the discrete pressure is defined on the primal grid, while the discrete velocity field and the density are defined on a face-based staggered dual grid. Then, the mass conservation equation, as well as the nonlinear convective terms in the momentum equation and the transport of kinetic energy in the energy equation are discretized explicitly, while the pressure terms appearing in the momentum and energy equation are discretized implicitly. Formal substitution of the discrete momentum equation into the total energy conservation equation yields a linear system for only one unknown, namely the scalar pressure. Here the equation of state is assumed linear with respect to the pressure. The enthalpy and the kinetic energy are taken explicitly and are then updated using a simple Picard procedure. Thanks to the use of a staggered grid, the final pressure system is a very sparse block five-point system for three dimensional problems and it is a block four-point system in the two dimensional case. Furthermore, for high order in space and piecewise constant polynomials in time, the system is observed to be symmetric and positive definite. This allows to use fast linear solvers such as the conjugate gradient (CG) method. In addition, all the volume and surface integrals needed by the scheme depend only on the geometry and the polynomial degree of the basis and test functions and can therefore be precomputed and stored in a preprocessing stage. This leads to significant savings in terms of computational effort for the time evolution part. In this way also the extension to a fully curved isoparametric approach becomes natural and affects only the preprocessing step. The viscous terms and the heat flux are also discretized making use of the staggered grid by defining the viscous stress tensor and the heat flux vector on the dual grid, which corresponds to the use of a lifting operator, but on the dual grid. The time step of our new numerical method is limited by a CFL condition based only on the fluid velocity and not on the sound speed. This makes the method particularly interesting for low Mach number flows. Finally, a very simple combination of artificial viscosity and the a posteriori MOOD technique allows to deal with shock waves and thus permits also to simulate high Mach number flows. We show computational results for a large set of two and three-dimensional benchmark problems, including both low and high Mach number flows and using polynomial approximation degrees up to p=4.

A High-Order Staggered Meshless Method for Elliptic Problems

We present a new meshless method for scalar diffusion equations, which is motivated by their compatible discretizations on primal-dual grids. Unlike the latter though, our approach is truly meshless because it only requires the graph of nearby neighbor connectivity of the discretization points ${\em x}_i$. This graph defines a local primal-dual grid complex with a virtual dual grid, in the sense that specification of the dual metric attributes is implicit in the method's construction. Our method combines a topological gradient operator on the local primal grid with a generalized moving least squares approximation of the divergence on the local dual grid. We show that the resulting approximation of the div-grad operator maintains polynomial reproduction to arbitrary orders and yields a meshless method, which attains $O(h^{m})$ convergence in both $L^2$- and $H^1$-norms, similar to mixed finite element methods. We demonstrate this convergence on curvilinear domains using manufactured solutions in two and th...

Reformulation and extension of the thrust network analysis

We address the Thrust Network Analysis (TNA), i.e. the methodology for modeling masonry vaults as a discrete network of forces in equilibrium with gravitational loads, first contributed by O’Dwyer and fully developed by Block and coworkers. Reducing the bias by the quoted authors in favor of a graphical interpretation of the method, we reformulate the original version of the TNA by discarding the dual grid and focusing only on the primal grid, thus significantly enhancing the computational performances. The proposed reformulation of the TNA is also extended by including horizontal forces in the analysis as well as holes or free edges in the vault. Furthermore, the coefficient matrices entering the solution scheme are obtained by assembling the separate contribution of each branch, thus avoiding the ad hoc node numbering and branch orientation required by alternative implementations. Numerical examples, some of which referred to vaults having a particularly complex geometry, show the effectiveness and robustness of the proposed approach in assessing the safety conditions of existing masonry vaults or in designing new ones.

New central and central discontinuous Galerkin schemes on overlapping cells of unstructured grids for solving ideal magnetohydrodynamic equations with globally divergence-free magnetic field

New schemes are developed on triangular grids for solving ideal magnetohydrodynamic equations while preserving globally divergence-free magnetic field. These schemes incorporate the constrained transport (CT) scheme of Evans and Hawley [34] with central schemes and central discontinuous Galerkin methods on overlapping cells which have no need for solving Riemann problems across cell edges where there are discontinuities of the numerical solution. These schemes are formally second-order accurate with major development on the reconstruction of globally divergence-free magnetic field on polygonal dual mesh. Moreover, the computational cost is reduced by solving the complete set of governing equations on the primal grid while only solving the magnetic induction equation on the polygonal dual mesh. Various numerical experiments are provided to validate the new schemes.

Geometrically defined basis functions for polyhedral elements with applications to computational electromagnetics

In the recent years, reformulating the mathematical description of physical laws in an algebraic form using tools from algebraic topology gained popularity in computational physics. Physical variables are defined as fluxes or circulations on oriented geometric elements of a pair of dual interlocked cell complexes, while physical laws are expressed in a metric-free fashion with incidence matrices. The metric and the material information are encoded in the discrete counterpart of the constitutive laws of materials, also referred to as constitutive or material matrices. The stability and consistency of the method is guaranteed by precise properties (symmetry, positive definiteness, consistency) that material matrices have to fulfill. The main advantage of this approach is that material matrices, even for arbitrary star-shaped polyhedral elements, can be geometrically defined, by simple closed-form expressions, in terms of the geometric elements of the primal and dual grids. That is why this original technique may be considered as a “Discrete Geometric Approach” (DGA) to computational physics. This paper first details the set of vector basis functions associated with the edges and faces of a polyhedral primal grid or of a dual grid. Then, it extends the construction of constitutive matrices for bianisotropic media.

A staggered semi-implicit discontinuous Galerkin method for the two dimensional incompressible Navier–Stokes equations

In this paper we propose a new spatially high order accurate semi-implicit discontinuous Galerkin (DG) method for the solution of the two dimensional incompressible Navier–Stokes equations on staggered unstructured curved meshes. While the discrete pressure is defined on the primal grid, the discrete velocity vector field is defined on an edge-based dual grid. The flexibility of high order DG methods on curved unstructured meshes allows to discretize even complex physical domains on rather coarse grids.Formal substitution of the discrete momentum equation into the discrete continuity equation yields one sparse linear equation system with four non-zero blocks per element for only one scalar unknown, namely the pressure. The method is computationally efficient, since the resulting system is not only very sparse but also symmetric and positive definite for appropriate boundary conditions. Furthermore, all the volume and surface integrals needed by the scheme presented in this paper depend only on the geometry and the polynomial degree of the basis and test functions and can therefore be precomputed and stored in a preprocessor stage, which leads to savings in terms of computational effort for the time evolution part. In this way also the extension to a fully curved isoparametric approach becomes natural and affects only the preprocessing step. The method is validated for polynomial degrees up to p=3 by solving some typical numerical test problems and comparing the numerical results with available analytical solutions or other numerical and experimental reference data.

A conservative adaptive wavelet method for the shallow‐water equations on staggered grids

AbstractThis article presents the first dynamically adaptive wavelet method for the shallow‐water equations (SWEs) on a staggered hexagonal C‐grid. Pressure is located at the centres of the primal grid (hexagons) and velocity is located at the edges of the dual grid (triangles). Distinct bi‐orthogonal second‐generation wavelet transforms are developed for the pressure and the velocity. These wavelet transforms are based on second‐order accurate interpolation and restriction operators. Together with compatible restriction operators for the mass flux and Bernoulli function, they ensure that mass is conserved and that there is no numerical generation of vorticity when solving the SWEs. Grid refinement relies on appropriate thresholding of the wavelet coefficients, allowing error control in both the quasi‐geostrophic and inertia–gravity wave regimes. The SWEs are discretized on the dynamically adapted multiscale grid using a mass and potential‐enstrophy‐conserving finite‐difference scheme. The conservation and error control properties of the method are verified by applying it to a propagating inertia–gravity wave packet and to rotating shallow‐water turbulence. Significant savings in the number of degrees of freedom are achieved even in the case of rotating shallow‐water turbulence. The numerical dissipation introduced by the grid adaptation is quantified. The method has been designed so it can be extended easily to the icosahedral subdivision of the sphere. This work provides important building blocks for the development of fully adaptive general circulation models.

A new Control Volume Finite Element Method for the stable and accurate solution of the drift–diffusion equations on general unstructured grids

We present a new Control Volume Finite Element Method with multi-dimensional Scharfetter–Gummel upwinding (CVFEM-SG) for the drift–diffusion equations. The method combines a conservative formulation of the carrier density continuity equations with an edge element lifting of the one-dimensional Scharfetter–Gummel edge currents into curl-conforming elemental currents. These elemental currents combine the upwind effect from all element edges and enable accurate computation of the flux on arbitrary surfaces inside the elements. In so doing, we obtain a formulation that is stable and accurate on general unstructured finite element grids. This approach sets our formulation apart from other methods, which require the control volumes to be topologically dual to the primal grid. Numerical studies of the CVFEM-SG for a suite of scalar advection–diffusion test problems confirm the accuracy and the robustness of the new formulation. Simulations of a PN diode and an n-channel MOSFET device demonstrate the performance of the method for the fully coupled drift–diffusion system.

Discrete Geometric Formulation of Admittance Boundary Conditions for Frequency Domain Problems Over Tetrahedral Dual Grids

In this work, it is shown how admittance boundary conditions for electromagnetic boundary value problems in the frequency domain can be formulated for the Discrete Geometric Approach. The details are presented for primal grids composed of tetrahedra and barycentric dual grids.

Constitutive Relations for Discrete Geometric Approach Over Hexahedral Grids

As a first step, we will investigate the extent of the consistency error for a hexahedral primal grid due to the use mass matrices of finite elements within the discrete geometric approach (DGA) as discrete counterparts of the constitutive relations. Such matrices are constructed by resorting to mixed elements edge and face vector functions. Then, we will propose a novel set of geometric vector basis functions for hexahedra that guarantee consistency and stability of the constitutive matrices for the DGA constructed according to the so-called energetic approach.

A new set of basis functions for the discrete geometric approach

By exploiting the geometric structure behind Maxwell’s equations, the so called discrete geometric approach allows to translate the physical laws of electromagnetism into discrete relations, involving circulations and fluxes associated with the geometric elements of a pair of interlocked grids: the primal grid and the dual grid. To form a finite dimensional system of equations, discrete counterparts of the constitutive relations must be introduced in addition. They are referred to as constitutive matrices which must comply with precise properties (symmetry, positive definiteness, consistency) in order to guarantee the stability and consistency of the overall finite dimensional system of equations. The aim of this work is to introduce a general and efficient set of vector functions associated with the edges and faces of a polyhedral primal grids or of a dual grid obtained from the barycentric subdivision of the boundary of the primal grid; these vector functions comply with precise specifications which allow to construct stable and consistent discrete constitutive equations for the discrete geometric approach in the framework of an energetic method.

Subgridding to Solving Magnetostatics Within Discrete Geometric Approach

We propose a recipe to construct a symmetric positive definite and consistent reluctance constitutive matrix to be used within discrete geometric approaches when the primal grid is generated by an enhanced subgridding of a generic hexahedral grid. We focus on a magnetostatic problem as working example.

Three-Dimensional Adaptive Central Schemes on Unstructured Staggered Grids

We present an explicit second-order finite volume generalization of the one-dimensional (1D) Nessyahu–Tadmor schemes for hyperbolic equations on adaptive unstructured tetrahedral grids. The nonoscillatory central difference scheme of Nessyahu and Tadmor, in which the resolution of the Riemann problem at the cell interfaces is bypassed thanks to the use of the staggered Lax–Friedrichs scheme, is extended here to a two-steps scheme. In order to reduce artificial viscosity, we start with an adaptively refined primal grid in three dimensions (3D), where the theoretical a posteriori result of the first-order scheme is used to derive appropriate refinement indicators. We apply those methods to solve Euler's equations. Numerical experimental tests on classical problems are obtained by our method and by the computational fluid dynamics software Fluent. These tests include results for the 3D Euler system (shock tube problem) and flow around an NACA0012 airfoil.

Base functions and discrete constitutive relations for staggered polyhedral grids

An electromagnetic problem can be discretized on a pair of interlocked primal–dual grids according to discrete geometric approaches like the Finite Integration Technique (FIT) or the Cell Method (CM). The critical aspect is however the construction of the discrete counterparts of the constitutive relations assuring stability and consistency of the overall discrete system of algebraic equations. Initially only orthogonal Cartesian grids where considered; more recently primal grids of tetrahedra and oblique prisms with triangular base can be handled. With this paper a novel set of edge and face vector functions for general polyhedral primal grids is presented, complying with precise specifications which allow to construct stable and consistent discrete constitutive equations in the framework of an energetic approach.

Explicit, Consistent, and Conditionally Stable Extension of FD-TD to Tetrahedral Grids by FIT

The finite difference-time domain algorithm is extended by means of finite integration technique to the case of tetrahedral primal grids. The resulting algorithm, unlike all previous attempts proposed in literature, is both explicit, consistent, and conditionally stable. It can be applied to both isotropic and anisotropic linear media. Numerical experiments demonstrate that the proposed algorithm leads to accurate approximations of the electromagnetic field.

3D adaptive central schemes: Part I. Algorithms for assembling the dual mesh

Central schemes are frequently used for incompressible and compressible flow calculations. The present paper is the first in a forthcoming series where a new approach to a 2nd order accurate Finite Volume scheme operating on Cartesian grids is discussed. Here we start with an adaptively refined Cartesian primal grid in 3D and present a construction technique for the staggered dual grid based on L ∞ -Voronoi cells. The local refinement constellation on the primal grid leads to a finite number of uniquely defined local patterns on a primal cell. Assembling adjacent local patterns forms the dual grid. All local patterns can be analysed in advance. Later, running the numerical scheme on staggered grids, all necessary geometric information can instantly be retrieved from lookup-tables. The new scheme is compared to established ones in terms of algorithmic complexity and computational effort.

Use of Barycentric Dual Grids for the Solution of Frequency Domain Problems by FIT

In this paper, the finite integration technique for the approximation of three-dimensional electromagnetic boundary value problems in the frequency domain is formulated for primal grids composed of either oblique parallelepipeds or oblique triangular prisms or tetrahedra, adopting as dual grid the barycentric subdivision of the primal one. Novel constitutive relations are in particular derived assuring solution matrices of symmetric and positive definite type.