Summary This paper describes a control-volume finite-element (CVFE) method incorporating linear triangular elements for the simulation of thermal multiphase flow in porous media. The technique adopts the usual finite-element shape functions to evaluate flow potentials at the control-volume boundaries and uses the conservation equations for each control volume. The main advantage of the CVFE method over the finite-difference method is in the representation of complex reservoir geometries. In addition, desirable features, such as local grid refinement for near-well resolution, can be achieved simply and consistently. The control-volume approach enforces local mass conservation and permits a direct physical interpretation of the resulting discrete equations. These are significant advantages over the classical Petrov-Galerkin or variational finite-element methods. The method was implemented in a general-purpose thermal simulator. Numerical examples compare the proposed method with five-point and nine-point finite-difference schemes in terms of grid-orientation effects and run time. The CVFE method was found to reduce grid-orientation effects significantly. At the same time, computational cost was much lower than for the nine-point scheme. The geometric flexibility of the method also is demonstrated. Introduction In many reservoir simulation problems, a flexible discretization method is extremely useful in the definition of complex reservoir geometries and discontinuities (such as faults) and in enhancing the resolution near the wells. The use of Cartesian grids with finite-difference methods has created difficulties and/or complexities in the definition of complex geometries or grid refinements.1-6 It is desirable to adopt the intrinsic grid flexibility of the finite-element method. However, combining upstream weighting with the usual finite-element method for the multiphase multidimensional flow problem presents difficulties. Although asymmetric weighting procedures like the Petrov-Galerkin method7 have been introduced to deal with the convective terms in the mixed convective-diffusive flow problems, such methods are in general not mass-conservative in the local sense. On the other hand, local mass conservation is a specific requirement of the control-volume methods. In addition, reservoir simulation problems can be very complex, involving multiphase mass and heat flow with interphase transfers and chemical reactions. The mass-conservative aspect of the control-volume methods is a distinct advantage in the programming and testing of these simulators. The CVFE method was proposed in computational fluid dynamics for solving the Navier-Stokes equations,8,9 where flexible gridding and local mass, momentum, and energy conservation are achieved. In this paper, a CVFE procedure for the reservoir flow equations is developed where flexible grid geometry is obtained without sacrificing the advantageous attributes of the control-volume finite-difference method. The derivation shows that the use of the perpendicular-bisection10 grid and the seven-point finite-difference method11 are special cases of this discretization method. Recently, Forsyth12 applied a CVFE method to the local-mesh-refinement problem by providing a smooth transition between the coarse and fine grids. As discussed in detail later, a proper choice of the triangular finite-element mesh is crucial to the reduction of grid-orientation effects. The construction of a CVFE grid by triangulation with one of the diagonals of each rectangle in a Cartesian grid (as in Ref. 12) will result in a five-point discretization scheme because the diagonal flow terms for this grid are identically equal to zero. The method results in a set of discretized conservative equations where the Jacobian construction for Newton's method and the upstream weighting of mobilities can be done in the usual way. For the incompressible single-phase flow problem, the method gives the same stiffness matrix as the Petrov-Galerkin weighted-residual finite-element method when linear shape functions are used. The criterion for maintaining positive transmissibility coefficients of a general anisotropic system also is derived. A number of examples are included to demonstrate the geometric flexibility, non-grid-orientation characteristics, and efficiency of the proposed method. p. 349-357

Giant reservoirs of the Middle East are crucial for the supply of oil and gas to the world market. Proper simulation of these giant reservoirs with long history and large amount of static and dynamic data requires efficient parallel simulation technologies, powerful visualization and data processing capabilities. This paper describes GigaPOWERS, a new parallel reservoir simulator capable of simulating hundreds of millions of cells to a billion cells with long production history in practical times. The new simulator uses unstructured grids. A distributed unstructured grid infrastructure has been developed for models using unstructured or complex structured grids. Unconventional wells such as maximum reservoir contact wells and fish-bone wells, as well as faults and fractures are handled by the new gridding system. A new parallel linear solver has been developed to solve the resulting linear system of equations. Load balancing issues are also discussed. A unified compositional formulation has been implemented. The simulator is designed to handle n-porosity systems. An optimization-based well management system has been developed by using mixed integer nonlinear programming. In addition to the core computational algorithms, the paper will present the pre- and post-processing software system to handle large amount of data. Visualization techniques for billions of cells are also presented. Introduction For many oil and gas reservoirs, especially large reservoirs in the Middle East, availability of vast amount of seismic, geological and dynamic reservoir data result in high-resolution geological models. But despite the many benefits of parallel simulation technology for large reservoirs, average cell size still remains in the order of hundreds of meters for large reservoirs. In order to fully utilize the seismic data, smaller grid blocks such as 25 to 50 meters in length are required. This size of grid blocks results in billion (Giga) cell models for giant reservoirs. In order to simulate such models with reasonable turnaround time, new innovations in the main components of the simulator such as linear equation solvers and equation of state computations are essential. Also, next generation pre- and post-processing tools are needed in order to build and analyze giga-cell models in practical times.

um .... ". A new switching criterion for adaptive-implicit reservoir simulation, based on the numerical stability of the local amplification matrix, is discussed. Adaptive-implicit methods can be used to achieve substantial savings in computing time and storage. 'The success of their application, however, depends critically on an appropriate switching criterion. 1be standard criterion discussed in the literature is the application of threshold changes in the primary variables, which are usually determined by some preliminary testing and chosen conservatively because of their empirical nature. Further, backward switching (implicit to explicit) is not normally feasible. With the switching criterion discussed in this paper, these disadvantages are eliminated. 'The method has been applied to a variety of black-oil reservoir simulation problems, including the first and second SPE Comparative Solution Projects and a water-injection problem. In all cases tested, the new criterion is shown to perform satisfactorily. For the waterflood problem in particular, the new criterion with backward switching was found to decrease computing time and storage significantly compared with the standard switching criterion. This criterion can easily be extended to other types of reservoir simulators. IntrocluotlonThe use of adaptive--implicit methods in reservoir simulation has been recognized as an optimal method of providing the best of explicit differencing methods and fully implicit methods. 'The former results in a substantially reduced system that is computationally inexpensive per timestep, but the numerical stability of such a system requires the use of small timesteps, which can be extremely restrictive. 'The latter does not have such a restriction on the timestep size but is much more expensive computationally per timestep be-cause all independent variables are treated implicitly, which results in a large system to be constructed and solved. Adaptive--implicit methods were pioneered by Thomas and Thurnau. I It was noted that for most reservoir simulation problems, a fully implicit treatment is usually not required and that if a fraction of the c0mputational cells are treated implicitly, sufficient computational stability is provided to allow for large timesteps of a size comparable to a fully implicit simulation. At that time, it was suggested that a threshold change of pressure and/or saturations could be used as the criterion for implicit switching of each individual cell . Further experience 2 showed that the use of threshold changes is sometimes not sufficient as a switching criterion. Thus, if threshold changes are used, conservative thresholds are required. For this reason, moreover, threshold changes cannot be used for switching implicit cells to. explicit cells and the optimal thresholds will differ from problem to problem and require adjustments or testing to determine their values. In view of this, an analytical criterion that is based on the numerical stability at a particular timestep is desirable.In this paper, an ...

A new parallel, black-oil-production reservoir simulator (Powers**) has been developed and fully integrated into the pre-and post-processing graphical environment. Its primary use is to simulate the giant oil and gas reservoirs of the Middle East using millions of cells. The new simulator has been created for parallelism and scalability, with the aim of making megacell simulation a day-to-day reservoir-management tool. Upon its completion, the parallel simulator was validated against published benchmark problems and other industrial simulators. Several giant oilreservoir studies have been conducted with million-cell descriptions. This paper presents the model formulation, parallel linear solver, parallel locally refined grids, and parallel well management. The benefits of using megacell simulation models are illustrated by a real field example used to confirm bypassed oil zones and obtain a history match in a short time period. With the new technology, preprocessing, construction, running, and postprocessing of megacell models is finally practical. A typical history-match run for a field with 30 to 50 years of production takes only a few hours.

Summary Simulation of fractured reservoirs with the dual-porosity/dual-permeabilityapproach involves discretization of the solution domain into two collocatedcontinua called the matrix and the fracture. The original idealized modelassumes that the matrix acts essentially as a source or sink to the fracture, which is the primary conduit for fluid flow. In multiphase flow situations, this idealization was found to be inadequate. Enhancements are needed torepresent the local matrix/fracture and matrix/matrix drainage and imbibitionprocesses. Attempts to represent these processes include the gravity-segregatedmodel, the subdomain model, the pseudofunction method, and thedual-permeability model. This work examines the mechanisms involved in gas/oilgravity drainage in terms of the block-to-block process. Current methods fortreating this problem are reviewed to identify deficiencies. A new approach isproposed in which these mechanisms can he represented properly in thefield-scale simulation of these reservoirs. The method involves the calculationof pseudo capillary potentials, which in an average sense (on a computationalblock basis) give the correct flow behaviors. These pseudos can be calculated apriori if a vertical equilibrium (VE) assumption can be made about the fluiddistribution in the matrix blocks. When the VE assumption is not valid, thepseudos can be determined from fine-grid simulations. Introduction The development of simulation tools that can capture the dominant flowprocesses and recovery mechanisms of naturally fractured reservoirs has beenthe subject of much attention and controversy over the last 15 years. Since theinception of the dual-porosity con-cept, many research papers on models thatuse the concept have appeared in the literature. The dual-porosity approachassumes that the fissured porous media can be represented by two overlappingcontinua called the fracture and the matrix. The fracture continuum consists ofthe interconnected network of fractures and/or so-lution vugs that constitutethe primary conduits for fluid flow. The matrix continuum consists of theintergranular pore space of the rock, which comprises the majority of thestorage in the pore space of the rock, which comprises the majority of thestorage in the reservoir. Early dual-porosity models include those of Kazemi etal. and Saidi. Saidi modeled a fractured reservoir by dividing it into sectorsin which the fracture was assumed to have infinite transmissi-bility. Thematrix was represented by several cylindrical matrix blocks that were suitablydistributed vertically and horizontally. These matrix blocks were gridded withboundary conditions imposed on them by the fracture, which was under thegravity-segregation assumption. Kazemi et al. discretized the fracturecon-tinuum into gridblocks and simulated fluid flow by a set of fracturemass-balance equations. The matrix was assumed to act as a source or sink tothe fracture, and the flow between the two continua was represented by a singlematrix/fracture transfer term. With the assumption that the matrix blocks wereisolated, the transfer term was constructed from a representative matrix blocklocated at the center of the gridblock. The total transfer rate for thegrid-block was thus the transfer rate for the representative matrix blockmultiplied by the total number of matrix blocks within the gridblock. Becausethe representative matrix blocks and the fracture were at the same depth, gravity effect on recovery from the matrix was not included in the transfercalculation. Most state-of-the-art dual-porosity simulators discretize thefracture continuum but with additional enhancements to handle the effects ofgravity on the transfer. Saidi et al. discussed the gas/oil gravity-drainage process in fracturedreservoirs in Iran. Because of the low viscous gradient in the fracture, liberated gas tends to percolate to form a secondary gas cap. As the gas zoneadvances, matrix blocks become surrounded by gas. The density differencebetween the gas in the fracture and the oil in the matrix becomes the maindriving force for oil recovery, which is limited by the capillary pressurebetween the two phases. The ultimate recovery from a single matrix block undergravity drainage depends on the balance between two forces: gravity, which is adirect function of the matrix block height and the density difference betweenthe two phases, and capillarity, which depends on the gas/oil interfacialtension (IFT). The two important phenomena associated with the gravity-drainage processesare reimbibition and matrix continuity. Reim-bibition refers to processes arereimbibition and matrix continuity. Reim-bibition refers to the re-entering ofdrained matrix oil into other matrix blocks, also known as the block-to-blockprocess. Matrix continuity describes the situation where the fractureseparating the matrix blocks is incomplete, resulting in some degree ofcontinuity of flow paths and potentials between them. These phenomena, as wewill see, have significant implications for the modeling phenomena, as we willsee, have significant implications for the modeling of gravity drainage withthe dual-porosity concept. Much of the more recent literature on dual-porosity modeling is devoted todeveloping enhancements for modeling the gravity effects in the transfercalculation. These can be roughly classified into four groups:gravity-segregated, subdomain, pseudofunction, and dual-permeability models. Ref. 15 discusses these models. Essentially, all the models try to account forgravity effects, with varying degrees of accuracy and complexity. None of themethods treats the reimbibition phenomenon, although the concept ofreimbibition is rather irrelevant to phenomenon, although the concept ofreimbibition is rather irrelevant to the dual-permeability model. All themodels except the dual-permeability model neglect capillary continuity. In thedual-permeability model, the matrix is assumed to be completely continuous. Thematrix oil recovery under gravity drainage reflects the gravity/capillarybal-ance for the entire matrix column and has nothing to do with the matrixblock height. This work examines the gravity-drainage process and associated phenomena, such as capillary continuity and reimbibition. A new method is then introducedthat can handle block-to-block interaction and is efficient for use in thefield-scale simulation of these reservoirs. The method involves the calculationof pseudo capillary potentials, which give the correct flow behaviors on acomputational block basis. Accuracy of the technique is verified by comparingits results with results from fine-grid simulations. Some comparisons withexisting models are also made.