Abstract
Reservoir simulation yields a system of linear algebraic equations, Ap=q, that may be solved by Richardson's iterative method, p(k+1)=p(k)+tkr(k), where r(k)=q-Ap(k) is the residual and t0, . . . tk are acceleration parameters. The incomplete factorization, Ka, of the strongly implicit procedure (SIP) yields an improvement of Richardson's method, p(k+1)=p(k)+tkKa−1r(k). Parameter a originates from SIP. The product of the L and U factors produced by SIP gives Ka=LU. The best values of the tk acceleration parameters may be computed dynamically by an efficient algorithm; the best value of a must be found by trial and error, which is not hard for only one value. The advantages of the method are (1) it always converges, (2) with the exception of the a parameter, parameters are computed dynamically, and (3) convergence is efficient for test problems characterized by heterogeneities and transmissibilities varying over 10 orders of magnitude. The test problems originate from field data and were suggested by industry personnel as particularly difficult. Dynamic computation of parameters is also a feature of the conjugate gradient method, but the iteration described here does not require A to be symmetric. Matrix Ka−1 A must be such that the real part of each eigenvalue is nonnegative, or the real part of each is nonpositive, but not both positive and negative. It is in this sense that the method always converges. This condition is satisfied by many simulator-generated matrices. The method also may be applied to matrices arising from the simulation of other processes, such as chemical flooding.
Introduction
The solution of a linear algebraic system, Ap=q, is a basic, costly step in the numerical simulation of a hydrocarbon reservoir. Many current solution methods are impractical for large linear systems arising from three-dimensional simulations or from reservoirs characterized by widely varying and discontinuous physical parameters. An iterative solution is described with these two main advantages:it is efficient for difficult problems andthe selection of iteration parameters is straightforward.
The method is Richardson's method applied to a preconditioned linear system. Matrix A may be symmetric or nonsymmetric. In the simulation of multiphase flow, it is usually nonsymmetric. Convergence behavior is shown for four examples. Two of these, Examples 3 and 4, were provided by an industry laboratory (Exxon Production Research Co.), and were suggested by personnel as especially difficult to solve; SIP failed to converge and only the diagonal method1 was effective. Convergence of Richardson's method is compared with the diagonal method using data from a laboratory run. The other two examples are: Example 1, a matrix not difficult to solve, generated from field data, and Example 2, a variant of a difficult matrix described by Stone.2 The easy matrix of Example 1 is included to show the performance of Richardson's method (with preconditioning) on a simple problem.
A Method for Imposing Surface Stress and Heat Flux Conditions in Finite-Difference Models with Steep Terrain