BACKGROUND.  The International Reservoir Simulation Research Institute was formed at Brigham Young University in July 1996. IRSRI is a new direction for a well-established research group. It is an attempt to make additional use of the resources of BYU's Advanced Combustion Engineering Research Center, ACERC, which was formed in 1986. One of ACERC's principal research areas has been the numerical simulation of combustion processes such as furnaces and turbines, a technology that has much in common with reservoir simulations. Since that inception IRSRI has gained the support of thirteen oil companies. The cost of these IRSRI memberships is $25,000 for the first year and $9,000 per year thereafter.

IRSRI's objective is to develop improved mathematical methods for petroleum reservoir simulation which will greatly enhance existing reservoir simulators. We hope to improve the speed and accuracy of traditional simulation methods to accommodate the large number of simulations required by the emerging technologies of geo-statistical reservoir descriptions, automatic field optimization, and automatic history matching. We also anticipate that these improvements will make full field reservoir simulation more practical on personal computers.

RESEARCH AREAS:  IRSRI's research includes projects which will improve the mathematics of reservoir simulation. The following projects are currently underway:

• Accurate Treatment of Complex Well Geometries
• Dynamic Gridding for Better Saturation Accuracy
• Flexible Cartesian Grids
• Prototype Simulator
• Curvilinear Orthogonal Grids

For further information on IRSRI, please e-mail irsri@byu.edu, call 801-378-3749, or write IRSRI, Brigham Young University, 45CTB, Provo, Utah 84602-4214.

Accurate Treatment of Complex Well Geometries.  A study was initiated at IRSRI's formation to investigate the use of a hybrid finite difference/analytical approach to simulation. An analytical solution based on well test theory was added to the traditional finite difference equations in an effort to eliminate the non-linearities which occur in the pressure solution as a result of the wells. The work was reported in "An Improved Method for Simulating Reservoir Pressures Through the Incorporation of Analytical Well Functions", SPE39065, presented at the Fifth Latin American and Caribbean Petroleum Engineering Conference and Exhibition held in Rio de Janeiro, Brazil, Sept. 1997.

The study showed that the hybrid approach had several distinct advantages:

• Improved Accuracy.  The elimination of near-singularities in the pressures near the wells from the finite difference equations did dramatically reduce their errors, as hypothesized. In seven, 2-D examples, the worst pressures errors were reduced by an average factor of 22.5 (2250%).
• Increased Speed.  This new approach eliminates the well flux values from the right-hand-side of the finite difference equations. Instead of a mostly sparse vector with occasional very large terms, at the well cells, the right-hand-side consists of relatively small, more consistent terms. It is anticipated that this new linear algebra will be much easier to solve. Indeed this was borne out by the relatively simple examples. Only about half the iterations were required for the new method as for the traditional approach.
• Elimination of Well Equations.  The analytical well functions actually predict the well bore pressures. There is no need for an empirical equation relating the well cell pressure to the well bore pressures which is used in the traditional approach. This advantage is of particular significance in situations where the traditional equations may not apply, e.g. complex well geometries or non-Cartesian grids.
• Flexible Well Placement.  The wells become independent of the finite difference equations. Therefore they need not be oriented with the grid. They can be vertical, horizontal, inclined, undulating, or any other shape that can be described as a series of straight-line segments. Wells need not be in the center of grid cells. They can be placed close together without regard for the number in intervening cells. Several well can even be in the same cell.
• Ease of Implementation.  The finite difference mathematics is much the same as in traditional simulators. Only the right-hand-sides of the linear algebraic equations are changed. Hence the method can be easily implemented in existing reservoir simulators.

These advantages of the new hybrid approach to the pressure solution are still believed to exist. However, recent multiphase results show that the principle errors in reservoir simulation result from numerical dispersion. Until the saturations (and concentrations) can be predicted more accurately, (See Dynamic Gridding for Better Saturation Accuracy) grossly improved pressure solutions are of little value. The principle value of this hybrid technology to existing reservoir simulators is its flexibility and accuracy in treating the wells.

Dynamic Gridding for Better Saturation Accuracy.  IRSRI's research into improved pressure solutions (See Accurate Treatment of Complex Well Geometries) has shown that the most significant reservoir simulation errors are the result of numerical dispersion of the saturations and concentrations. Several approaches have been investigated by IRSRI without great success. Numerical dispersion has been the bane of numerical methods for solving partial differential equations since their beginning. Much work has been done within the Petroleum Industry and elsewhere to overcome the problem. The most viable solution to date seems to be simply the reduction of grid block sizes. Although smaller grid blocks do reduce numerical dispersion, they also result in increased computation times. To minimize the increased number of blocks, and their inherent cost, many investigators have considered a dynamic grid in which the grid is refined only in the areas where dispersion is severe. Grid blocks are constantly being introduced into the simulation and removed as the dispersion locations change. However, this type of dynamic gridding has not been widely accepted, because the substantial additional computation required to generate the constantly changing grids.

IRSRI's most recent attempt to increase saturation accuracy is through a dynamic gridding scheme which requires no extra computation. This is accomplished by solving the traditional finite difference equations, but instead of solving for saturation as the unknown, the grid locations are treated as unknowns. That is, rather than finding the saturations at particular grid points, the location of grid points are determined at particular saturations. Hence the grid becomes a series of constant saturation lines which move through the reservoir as the saturations change. Grids are automatically refined in areas where saturation gradients are large, and substantial numerical dispersion would be anticipated.

In one-dimension, the scheme is simple and works well. In the classical Buckley-Leverett problem, grid points cluster around the location of the flood front, resulting in dramatically sharper injection profiles. In multiple dimensions, however, the method is complicated by the fact that a particular saturation exists at many points throughout the reservoir. There are additional degrees of freedom available in defining the dynamic grid. One gridding scheme which seems to work well in multiple dimensions is to use the constant potential lines and flow lines to constitute the grid. Such a grid is orthogonal so that the original equations still apply. Moreover such a grid reduces to a series of one-dimensional flow problems. A two dimensional grid block is comprised of two sides which are flow lines and two sides which are constant potential lines. There is no flow across the two flow lines; flow occurs only between the constant potential sides. A three-dimensional grid block has four no-flow sides, and again the flow is one-dimensional between the constant potential sides. The new dynamic grid requires only solving a series of one-dimensional Buckley-Leverett-like problems down a series of stream tubes. The method is very fast and substantially reduces numerical dispersion.

The use of these dynamic grids comprised of stream tubes and constant potential surfaces to solve the saturation equations results in an algorithm which is very similar to characteristic methods, front tracking methods, and stream tube methods, which have been popularized recently. However these other techniques require the tracking of two types of surfaces, both saturation contours and flood fronts. The new method calculates only the saturation contours, the fronts take care of themselves. This can be a considerable advantage in reservoirs of complex heterogeneity where many flood fronts may form and disappear throughout the simulation. The new method can also include all physical phenomena for which finite difference equations can be written. The other methods have difficulties with some phenomena, e.g. capillary pressure and gravitational effects.

Flexible Cartesian Grids.  Cartesian coordinates have been widely used over the years for reservoir simulation because of their simplicity, speed, and accuracy. However, Cartesian grids lack flexibility in representing complex reservoir geometries, including both external boundaries and internal faults. As a result, much work has been done in recent years to develop more flexible gridding methods. The most promising of these are unstructured grids, which result in more complex linear algebra and increased computation time.

IRSRI has developed, and is assessing, mathematical boundary conditions for Cartesian coordinates which represent complex boundary conditions more accurately. This work will reduce the need for other, more complex grid types. The work provides Cartesian finite difference formulations which accurately represent boundaries passing obliquely through the finite difference cells. The method involves:

• The inclusion of all cells which contain a part of the reservoir. Hence, the new method usually includes more cells than traditional "stair-step" Cartesian boundaries.
• Modification of cell porosities to accurately reflect reservoir pore volumes.
• Modification of linking permeabilities to accurately reflect cell interface areas contained in the reservoir.

Internal faults which do not coincide with block boundaries are divided in two, each representing the reservoir on different sides of the fault. This alters the structure of the linear algebra matrices. Hence accurate treatment of complex external boundaries can be obtained using any existing Cartesian finite difference reservoir simulator. No modification is required. However, the treatment of internal faults will require changes to the simulator unless it features non-neighbor connection capabilities.

Prototype Simulator.  Much of IRSRI's research is sufficiently mature that it can be incorporated into existing reservoir simulators and into actual field simulation studies. To demonstrate their utility, a prototype simulator is being written. It will include:

• Accurate Treatment of Complex Well Geometries
• Dynamic Gridding for Better Saturation Accuracy
• Flexible Cartesian Grids

It is hoped that this simulator will start a new generation in reservoir simulation. Its speed and accuracy should accommodate very large problems. Its accuracy should accommodate very complex problems with modest grid sizes.

Curvilinear Orthogonal Grids.  Because of the inflexibility of Cartesian reservoirs in representing complex reservoir geometries, much work has been done in recent years to develop more flexible gridding systems. Corner Point grids have been widely used. PEBI grids have been touted. Other unstructured grids have been proposed. However, each of the new grid systems suffers from either reduced accuracy or reduced computation speeds. In each instance, there appears to be a price associated with the new grid's flexibility.

Curvilinear Orthogonal Grids provide a unique exception. Curvilinear orthogonal grids result from a conformal mapping of a complex reservoir geometry into a square (2-D) or a cube (3-D). Simulations are done using a Cartesian grid in the square or cube. The results are then transformed back into the curvilinear coordinates of the physical system for viewing and analysis. If the conformal mapping is exact, no accuracy is lost in the simulation. Since the simulations are actually done in Cartesian coordinates, no speed is lost either. In fact simulations are often slightly faster because inactive cells are eliminated.

The advantages of curvilinear orthogonal grids have long been recognized. However, there has not been a fast and robust method for conformally mapping arbitrary geometries into a rectangle. The principle breakthrough of this work has been the development of such an algorithm¹. It is now incorporated into a user-friendly, stand-alone gridding system, COG. The 2-D, curvilinear orthogonal grids generated by COG can be used in 3-D simulations by using the same grid in each layer, in much the same way that Cartesian grids are used. COG is a preprocessor program for reservoir simulation data. It generates a curvilinear orthogonal grid and calculates the necessary permeability and porosity multipliers that must be added to the physical data to accomplish the change in coordinates. It can also be used to interpolate the grid property values, such as permeabilities and porosities, which correspond the new grid. These data then can be used with any simulator capable of running simulations with Cartesian coordinates.

Work is currently underway to develop a similar computer program to generate 3-D curvilinear orthogonal grids.