Comparison of two finite‐difference approaches and five solution algorithms for 3‐D resistivity modeling

1994 ◽  
Author(s):  
Burkhard Wurmstich ◽  
Klaus Spitzer
1985 ◽  
Vol 25 (04) ◽  
pp. 535-542 ◽  
Author(s):  
A. Behie ◽  
D. Collins ◽  
P.A. Forsyth ◽  
P.H. Sammon

Abstract A fully coupled treatment of oil wells that are completed in more than one zone results in a bordered matrix. This paper develops solution algorithms that incorporate paper develops solution algorithms that incorporate existing direct and iterative (incomplete LU) solutions in a straightforward manner. Timings in scalar and vector modes on the Cray for a typical reservoir simulation problem are presented. problem are presented. Introduction Numerical simulation of oil reservoirs requires the solution of coupled sets of highly nonlinear partial differential equations. These equations represent the conservation of oil, gas, water, and energy. It usually is necessary to solve from 3 to 10 coupled equations per finite-difference cell. The equations usually are discretized by use of a nearest-neighbor coupling in space and a fully implicit timestep scheme. The resulting set of nonlinear algebraic equations then is solved by Newtonian iteration., Clearly, simulation of large systems requires effective solution of the Jacobian matrix. Many practical reservoir simulation problems involve multiblock wells or fractures. These situations arise when a well is completed in several layers, and consequently the wellbore penetrates several finite-difference cells. Each conservation equation in a cell penetrated by a well will have a source term of the form .....................................(1) where qjt is the mass influx of component k (resulting from the well), Xk is the mobility of component k, 1 pi is the pressure in cell i, and pi, is the unknown wellbore pressure in well j. pressure in well j. To specify the wellbore pressure, pi, an additional equation is required. This extra equation as generally a constraint op the total flow into the well - This constraint is of the form .....................................(2) where qJt. is the total specified fluid flow into well j, Nc, is the total number of components, and is the set of cell numbers penetrated by well j. Because several cells are connected to the same well, there is now an extra degree of coupling between these cells through the well-bore pressure. This coupling generally will not be consistent with the coupling produced by the usual finite-difference molecule. If the well pressures, pjw, are treated explicitly, or are lagged one iteration, convergence difficulties or stability limitations often result. 7 Fully coupled treatment of multiblock wells gives rise to a bordered matrix. We develop various methods to solve these systems. These methods are specifically designed for the block-banded systems arising from fully implicit thermal problems, although similar methods can be used for single-component systems The iterative methods are extensions of the incomplete factorization techniques (ILU), and a direct method is presented for comparison. Existing solution routines can be modified easily to solve the bordered system. Solution of the Bordered Matrix The standard approach to solving fully implicit, fully coupled multiblock wells (or fractures) is to order the unknowns so that those connected with flow in the reservoir (cell pressures, saturations, etc.) appear first in the solution vector. The unknowns connected with the well (well pressures) are placed last in the solution vector. This produces a bordered Jacobian matrix (see Fig. 1). The upper left portion of the matrix has the usual incidence matrix for the Jacobian of nearest-neighbor finite-difference discretization. The incidence matrix for the Jacobian is a matrix with entries zero if the Jacobian elements are zero, and with entries one if the Jacobian elements are nonzero. The border of columns on the upper right of Fig. 1 contains derivatives of the source terms (Eq. 1) with respect to the wellbore pressure The border of rows on the lower left contains derivatives of the constraint equations (Eq. 2) with respect to reservoir variables (i.e., cell pressures). The block on the lower right contains derivatives of the constraint equations with respect to the wellbore pressures and is diagonal. The number of extra columns and rows is proportional to the number of fully coupled wells (or fractures). Although the incidence matrix of the reservoir flow portion of the matrix is symmetric, the incidence matrix of portion of the matrix is symmetric, the incidence matrix of the borders is not necessarily symmetric. George discusses three possible block factorizations of sparse, linear systems. The algorithm used here is based on his second factorization. SPEJ P. 535


Geophysics ◽  
1978 ◽  
Vol 43 (5) ◽  
pp. 930-942 ◽  
Author(s):  
Irshad R. Mufti

Highly efficient finite‐difference resistivity modeling algorithms which yield accurate results are put forward. The given medium is discretized and divided into rectangular blocks by using a very coarse system of vertical and horizontal grid lines, whose distance from the source(s) increases logarithmically. Expressions are derived to compute the longitudinal conductance and transverse resistance associated with each of these blocks for a parallel‐layer medium followed by a generalized treatment to accommodate arbitrarily shaped structures. Since the values of Dar Zarrouk parameters are derived from the exact resistivity distribution of the given medium, fine features such as a thin but anomalously resistive bed which ordinarily would be missed entirely in coarse discretization can be taken into account. Further reduction in the size of the model is achieved by making use of a symmetry wherever possible. In most cases the computation of the potential field which involves the inversion of a small sparse matrix requires about 0.5 sec of computer time. Moreover, changes in geology affect neither the size nor the zero structure of the matrix. Therefore, when more than one model is to be computed, the factorization of the matrix can be done symbolically only once for all models, followed by numeric factorization for each individual model. The coarse grid algorithm was applied to a number of horizontally layered models involving a point source. The results obtained for each model were in excellent agreement with the corresponding analytical data. Finite‐difference investigation of the potential field for two‐dimensional structures and a line source dipole indicates that as long as one is interested only in the evaluation of the Schlumberger‐type apparent resistivity curves, the line‐source results may be a much better approximation to the corresponding point‐source data than is commonly believed.


Geophysics ◽  
2008 ◽  
Vol 73 (3) ◽  
pp. F135-F142 ◽  
Author(s):  
Erhan Erdoğan ◽  
Ismail Demirci ◽  
Mehmet Emin Candansayar

We incorporate topography into the 2D resistivity forward solution by using the finite-difference (FD) and finite-element (FE) numerical-solution methods. To achieve this, we develop a new algorithm that solves Poisson’s equation using the FE and FD approaches. We simulate topographic effects in the modeling algorithm using three FE approaches and two alternative FD approaches in which the air portion of the mesh is represented by very resistive cells. In both methods, we use rectangular and triangular discretization. Furthermore, we account for topographic effects by distorting the FE mesh with respect to the topography. We compare all methods for accuracy and calculation time on models with varying surface geometry and resistivity distributions. Comparisons show that model responses are similar when high-resistivity values are assigned to the top half of the rectangular cells at the air/earth boundary with the FE and FD methods and when the FE mesh is distorted. This result supports the idea that topographic effects can be incorporated into the forward solution by using the FD method; in some cases, this method also shortens calculation times. Additionally, this study shows that an FD solution with triangular discretization can be used successfully to calculate 2D DC-resistivity forward solutions.


Geophysics ◽  
1996 ◽  
Vol 61 (5) ◽  
pp. 1301-1307 ◽  
Author(s):  
Shengkai Zhao ◽  
Matthew J. Yedlin

Two basic refinements of the finite‐difference method for 3-D dc resistivity modeling are presented. The first is a more accurate formula for the source singularity removal. The second is the analytic computation of the source terms that arise from the decomposition of the potential into the primary potential because of the source current and the secondary potential caused by changes in the electrical conductivity. Three examples are presented: a simple two‐layered model, a vertical contact, and a buried sphere. Both accurate and approximate Dirichlet boundary conditions are used to compute the secondary potential. Numerical results show that for all three models, the average percentage error of the apparent resistivity obtained by the modified finite‐difference method with accurate boundary conditions is less than 0.5%. For the vertical contact and the buried sphere models, the error caused by the approximate boundary condition is less than 0.01%.


Geophysics ◽  
1996 ◽  
Vol 61 (6) ◽  
pp. 1616-1623 ◽  
Author(s):  
Shengkai Zhao ◽  
Matthew J. Yedlin

We use the multidomain Chebyshev spectral method to solve the 3-D forward direct current (dc) resistivity problem. We divided the whole domain into a number of subdomains and approximate the potential function by a separate set of Chebyshev polynomials in each subdomain. At an interface point, we require that both the potential and the flux be continuous. Numerical results show that for the same accuracy the multidomain Chebyshev spectral method is 2 to 260 times faster than the finite‐difference method.


Sign in / Sign up

Export Citation Format

Share Document