Gridding in reservoir simulation

The aim of gridding in reservoir simulation is to turn the geological model of the field into a discrete system on which the fluid flow equations can be solved.

Overview

The basic structure of an oil reservoir is a set of geological horizons representing bedding planes. The reservoir may contain faults, at which the strata are displaced. It is usually possible to identify many more layers in the geological model than it is practical to include in reservoir flow simulation, so some upscaling of rock properties will normally be carried out. Even after this process, the geology to be represented is rarely homogeneous at the scale of the simulation grid.

Two related issues are involved in choosing a grid for reservoir simulation:

• Accuracy with which the geological description of the reservoir is matched
• Discretization of the flow equations

In a classical finite-difference scheme, the point values of pressures and saturations are used as solution variables, and the differential operators that appear in the fluid-flow equations may be expanded as difference expressions of these point values to some order. An alternative approach is to use an integral finite-difference or finite-volume method in which the fluid-flow equations are integrated over a set of cell volumes. This yields a set of equations in which the mass conservation conditions for the fluid in the simulation cell volumes are related to the flows through the interfaces between those cell volumes. Rock properties such as porosity are assumed constant over the cell or controlled volume. This yields a discretization scheme which is conservative (each outflow from one cell is an inflow to another) and for which the fluid in place may be obtained straightforwardly. The mass conservation equations for a timestep from T to T + ΔT then become:

....................(1)

where

• V_pa is the pore volume of cell a
• m_ca is the density of conserved component c in cell a
• Q_ca is the injection or production rate of component c because of wells
• F_cpab is the flow rate of component c in phase p from cell a to its neighbor b

In general, the flows F_cpab may involve the solution values of a number of cells, the number of cells involved defining the stencil of the numerical scheme. The linear pressure dependence of flows given by Darcy’s law leads to an expression of the type:

....................(2)

M_cpax is the mobility of component c in phase p for the contribution to the flow between a and x, given by x_cp.K_rp/μ_p, where

• x_cp is the concentration of component c in phase p
• K_rp is the relative permeability of phase p
• μ_p is the viscosity of phase p

This is often set to an upstream value of the mobility, depending upon the sign of the potential difference.

ΔФ_pax is the potential difference of phase p between cell a and cell x, which includes pressure, gravity and capillary pressure contributions:

....................(3)

The constant coefficients of mobility and potential difference products, T_ax , are commonly termed the transmissibilities.

When the flows between two cells a and b can be expressed as a function of the solution values in just those two cells, so that the summation over cells includes just x = b, the flow expression takes a two-point form. The flow expression then takes a simple form:

....................(4)

When solution values from other cells are required, the flow takes a multipoint form.

Other options for discretization are available, such as Galerkin finite elements and mixed finite-element. It is sometimes possible to cast a finite-element Galerkin discretization into the upstreamed transmissibility-based form.

Regular cartesian grids

A simple 3D grid is the regular Cartesian grid (Fig. 1). Cells in such a grid may be simply identified using their (i,j,k) index values.

• 

  Fig. 1 – Regular Cartesian grids.

Each of the grid elements will be assigned a single permeability or porosity value. In this case, it is possible to obtain the transmissibility value as a harmonic average:

....................(5)

where

• cell b is the neighbor to cell a in some direction
• K is the cell permeability in that direction
• A is the area of the cell orthogonal to the direction of flow
• d the dimension of the cell in that direction

Such a two-point transmissibility assumes a permeability tensor with primary axes aligned along the grid axes.

Although regular grids are normally defined in normal Cartesian coordinates, it is also possible to use an (r, Ф, z) radial system. The resulting grid is cylindrical and is important for the special case of near-well studies dominated by radial inflow. For a 3D system, regular grids yield seven-point schemes, in which the flow equations for a cell involve solution values for just the cell and its six neighbors. Not all the elements in the grid need represent active solution variables in the simulation. Some cells may be inactive, representing volumes of the reservoir with zero porosity. Such inactive cells are usually compressed out of reservoir simulation solution arrays prior to the memory and time-intensive flow solution stage, and enable reservoirs with irregular boundaries to be represented within extended simulation grids.

The horizons that delimit rock strata are generally not horizontal, but are dipped, curved, or faulted. Unless extremely fine, a true regular grid that is orthogonal in all three axes will be unable to assign rock properties accurately to cell volumes. Such a layer-cake structure can be used, but will generally misalign property values (Fig. 2) in which the orthogonal grid provides a rather poor match to the dipping strata represented by the shaded layers. However, it is possible that improving computer power will bring such rasterized grids to a level of refinement at which a sufficiently good representation may be obtained.

• 

  Fig. 2 – Orthogonal grid used to represent dip.

Dip-normal geometry

A simple variation of a regular grid, in which the regular grid is rotated to bring the layers of cells into alignment with the bedding planes. Such a description would only suit a reservoir with a single, constant angle of dip. As geological descriptions have improved, fewer and fewer model reservoirs are found to fit this simple pattern, and something more flexible is required.

Block-center geometry

A simple model in which transmissibility between blocks is calculated on the basis of linear interpolation between the center values of the cells. This is a simple way of representing variable dip, but is difficult to represent graphically in a consistent way. Pore volumes are calculated on the basis of a series of flat regular cells with variable depths (Fig. 3a), but transmissibilities are calculated on the basis of interpolated values (Fig. 3b). The areal grid is rectangular.

• 

  Fig. 3 – (a) Obtaining pore volumes; (b) obtaining transmissibilities.

Thus, for the pair of cells illustrated,

....................(6)

where A is the average area over which flow occurs and c is a dip correction given by cos²θ, where θ is the angle of dip of a line joining the cell centers to the horizontal. Such a block-center option is suitable for unfaulted reservoirs and is commonly supplied as a simulator option.

Hexahedral grids

Further improvements in geological modeling threw an emphasis on describing faults, and made it important to distinguish depth displacements due to dip and faulting. This is difficult in block centre geometry in which the cell is positioned by its centre depth and Δx, Δy, Δz dimensions. To define faulting more precisely it is useful to define the position of grid cell by its corner point locations. A hexahedral shape with eight corners and bilinear planes as surfaces then describes the cell geometry. Faults, both vertical and inclined, may be described precisely (Fig. 4). Such grids are often called corner-point grids.

• 

  Fig. 4 – Hexahedral grid system.

In both the dipped and general hexahedral grids, the orthogonality of a completely rectangular grid no longer exists, and the result is that the two-point property of the flows between the cells is lost—the flow between cell a and cell b is not just a function of the solution values of cells and a and b. Typically, the result is a 27-point scheme in three dimensions. However, if the grid distortion is mild, it may be possible to ignore some additional couplings and use a low-order transmissibility scheme. This is normally done for extra couplings introduced by dip angles, which are often small.

Although this corner-point description handles the fault issue, the basic coordinate system remains a regular grid (i.e., the grid is structured). Fitting such a basically regular system to the irregular shapes of a reservoir remains a difficulty that may be solved in two basic ways:

• By distorting the grid and fitting the cells into the geometry

or

• By truncating the grid to the reservoir position

Multiple-domain hexahedral grids

In some cases, a single structured grid system cannot match the overall structure of a reservoir, so a block grid or domain-based grid is used. This consists of a number of subgrids, each with a local regular (i,j,k) structure, but linked together to model the entire reservoir. The block hexahedral system gives rise to multiple (i,j,k) indexing systems—(i,j,k,l), where the l index specifies the grid system. These comprise a series of regular grids. Such regular gridding systems have advantages for upscaling and downscaling—for example, a natural coarsening of a regular grid may be simply defined by grouping sets of coordinates in each direction.

Grid refinement

A common requirement in reservoir simulation is an increased level of detail around an item of interest such as a well. This is frequently obtained in structured grids by local grid refinement, replacing a set of cells in the original grid by a finer grid (Fig. 5). The inserted grid may be Cartesian (center) or radial (upper left). Local refinement may be regarded as a form of multiple domain structured grid, in that it consists of a number of linked structured grids. Flows at the edges of local refinements generally take a multipoint form.

• 

  Fig. 5 – Local grid refinement.

Unstructured grids

The problems involved in using a regular structured grid to represent reservoir geometry can be avoided by using an unstructured grid. This is constructed around a set of solution points that need have no particular indexing scheme. These points may be triangulated into a mesh of triangles or tetrahedrons. A control volume is constructed around the nodes of the resulting mesh to define the simulator cell volumes. The perpendicular bisector (PEBI) method introduced into reservoir simulation by Heinemann used a technique also known as a Voronoi grid. Starting from any set of solution points, the PEBI cell volumes are defined by the perpendicular bisection planes between these points. The resulting control volume is defined by the perpendicular planes—it is the set of points closer to the node than any other. This is shown in Fig 6, in which the bisectors to the heavy lines joining the solution points enclose the control volume, represented by the shaded area. The grid is locally orthogonal, and the desirable property of two point flows is obtained. The actual cell volumes may have a variety of shapes, depending on the exact placement of the solution points, but are typically hexagonal in two dimensions. Grid refinement occurs naturally in areas where solution points are closely spaced.

• 

  Fig. 6 – Unstructured grids may result in hexagonal cells.

The two-point property is not naturally preserved in anisotropic reservoirs, although it can be regained by transforming to a K-orthogonal grid in which the geometry is transformed so that K_n is parallel to the vector joining the solution nodes, where K is the permeability tensor and n is the normal to the cell volume surface. For nonisotropic cases in which the grid is not K-orthogonal, the flows will be functions of the solution values in more than two cells, as in the general hexahedral case.

An unstructured grid may be defined in two dimensions, and then applied to each layer of a reservoir model, so that a typical cell is a hexagonal prism. This is sometimes termed a 2½D unstructured grid. Alternatively, a full 3D unstructured grid may be defined. The 3D approach is most effective when applied to model a local structure such as a branched well.

Unstructured grids yield an elegant and flexible grid description. However, the ability to identify cells by a simple set of indices is lost, and items such as wells need to be positioned in true space terms. The systems of linear equations generated by unstructured grids are also commonly regarded as more difficult to solve than those produced by structured grids. However, it may be more true to say that optimal solution schemes are simpler to find for structured grids, where the row and plane order provides a natural ab initio solution variable-ordering scheme.

Truncated regular grids

The truncated approach fits in well with the rectangular grids used in geological modeling. A simple rectangular grid is always used in the areal direction, but faults may subdivide the rock volume in a given column. The areal grid is not modified to match the faults. Thus the two marked volumes in Fig. 7 represent different cells, but may have the same i, j indices, so this creates a multiple-domain grid. A disadvantage is the more complex shape of cells at the edge of the grid. Transmissibilities for such cells may need to be calculated numerically. Apart from the truncated cells, all the grid cells are hexahedra that are rectangular in plan.

• 

  Fig. 7 – A fault creating two cell volumes in a truncated grid.

Other gridding systems

Triangular or tetrahedral grids

The underlying solution points of a PEBI mesh can be linked together into a Delaunay triangulation. In 2D, this creates triangles, and in 3D it creates tetrahedra. One option would clearly be to use triangular or tetrahedral cells directly and associate cell volumes with these. This technique is rather rarely used in reservoir simulation. Partly this may be historical, but the Delaunay triangulation is rather less stable under grid changes than a Voronoi grid, and triangulation can more often lead to "sliver" cells with a high surface area but a small volume.

Curvilinear grid systems

In some special cases, a transformed coordinate system may be used, based around an expected flow pattern. Such grids are not well adapted to represent geological data, and have been used less frequently as more detailed reservoir descriptions have become available.

Future directions

Two themes emerge from current trends in reservoir simulation gridding. The increasing sophistication of data preparation and solver technology indicates a move towards unstructured grids as a general method of solving the flow equations for a given reservoir simulation problem. On the other hand, reservoir simulation is increasingly seen as part of a decision-making process rather than an isolated activity, so the ability to map easily onto the generally regular data structures used in seismic and geological modeling becomes an important issue. In this role, structured grids may have advantages of simplicity and scalability.

An ideal is to separate the construction of the flow-simulation grid from the description of the reservoir geometry. This ties in with a further ideal, inherent in many discretization schemes, that the scale of the simulation grid should be below the scale of the problem structure.

For more complex shape-dominated problems, the unstructured approach looks general and flexible, providing that the data-handling and cell-identification methods can be moved to true x,y,z space preprocessing software.

Nomenclature

d	=	dimension in direction of flow
K	=	permeability, md
ρ	=	density, mass/volume or molar density, mols/volume unless noted otherwise

Noteworthy papers in OnePetro

Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read

External links

Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro.

See also

Reservoir simulation

Upscaling of grid properties in reservoir simulation

Reservoir simulation applications

Geostatistical reservoir modeling

Geomechanics in reservoir simulation

Phase behavior in reservoir simulation

PEH:Reservoir_Simulation