Title of Invention  "MULTISCALE FINITEVOLUME METHOD FOR USE IN SUBSURFACE FLOW SIMULATION" 

Abstract  A multiscale finitevolume (MSFV) method to solve elliptic problems with a plurality of spatial scales arising from single or multiphase flows in porous media is provided. The method efficiently captures the effects of small scales on a coarse grid, is conservative, and treats tensor permeabilities correctly. The underlying idea is to construct transmissibilities that capture the local properties of a differential operator. This leads to a multipoint discretization scheme for a finitevolume solution algorithm. Transmissibilities for the MSFV method are preferably constructed only once as a preprocessing step and can be computed locally. 
Full Text  MULTISCALE FINITEVOLUME METHOD FOR USE IN SUBSURFACE FLOW SIMULATION RELATED APPLICATIONS This application is a continuationinpart of copending U.S. Patent Application Serial No. 10/383,908, entitled "MultiScale FiniteVolume Method for Use in Subsurface Flow Simulation", filed on March 6, 2003 and is a continuationinpart of copending application entitled "MultiScale FiniteVolume Method for Use in Subsurface Flow Simulation", filed on September 22, 2004, which is a continuation of copending U.S. Patent Application Serial No. 10/383,908, entitled "MultiScale FiniteVolume Method for Use in Subsurface Flow Simulation", filed on March 6, 2003. TECHNICAL FIELD The present invention relates generally to subsurface reservoir simulators, and more particularly, to those simulators which use multiscale physics to simulate flow in an underground reservoir. BACKGROUND OF THE INVENTION The level of detail available in reservoir description often exceeds the computational capability of existing reservoir simulators. This resolution gap is usually tackled by upscaling the finescale description to sizes that can be treated by a fullfeatured simulator. In upscaling, the original model is coarsened using a computationally inexpensive process. In flowbased methods, the process is based on singlephase flow. A simulation study is then performed using the coarsened model. Upscaling methods such as these have proven to be quite successful. However, it is not possible to have a prior estimate of the errors that are present when complex flow processes are investigated using coarse models constructed via these simplified settings. Various fundamentally different multiscale approaches for flow in porous media have been proposed to accommodate the finescale description directly. As opposed to upscaling, the multiscale approach targets the full problem with the original resolution. The upscaling methodology is typically based on resolving the length and timescales of interest by maximizing local operations. Arbogast et al. (T. Arbogast, Numerical subgrid upscaling of two phase flow in porous media, Technical report, Texas Institute for Computational and Applied Mathematics, The University of Texas at Austin, 1999, and T. Arbogast and S.L. Bryant, Numerical subgrid upscaling for waterflood simulations, SPE 66375, 2001) presented a mixed finiteelement method where finescale effects are localized by a boundary condition assumption at the coarse element boundaries. Then the smallscale influence is coupled with the coarsescale effects by numerical Greens functions. Hou and Wu (T. Hou and X.H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comp. Phys., 134:169189,1997) employed a finiteelement approach and constructed specific basis functions which capture the small scales. Again, localization is achieved by boundary condition assumptions for the coarse elements. To reduce the effects of these boundary conditions, an oversampling technique can be applied. Chen and Hou (Z. Chen and T.Y. Hou, A mixed finite element method for elliptic problems with rapidly oscillating coefficients, Math. Comput., June 2002) utilized these ideas in combination with a mixed finiteelement approach. Another approach by Beckie et al. (R. Beckie, A.A. Aldama, and E.F. Wood, Modeling the largescale dynamics of saturated groundwater flow using spatial filtering, Water Resources Research, 32:12691280,1996) is based on large eddy simulation (LES) techniques which are commonly used for turbulence modeling. Lee et al. (S.H. Lee, L.J. Durlofsky, M.F. Lough, and W.H. Chen, Finite difference simulation of geologically complex reservoirs with tensor permeabilities, SPERE&E, pages 567574,1998) developed a fluxcontinuous finitedifference (FCFD) scheme for 2D models. Lee et al. further developed a method to address 3D models (S.H. Lee, H. Tchelepi, P. Jenny and L. Dechant, Implementation of a flux continuous finitedifference method for stratigraphic, hexahedron grids, SPE Journal, September, pages 269277, 2002). Jenny et al. (P. Jenny, C. Wolfsteiner, S.H. Lee and L.J. Durlofsky, Modeling flow in geometrically complex reservoirs using hexahedral multiblock grids, SPE Journal, June, pages 149157, 2002) later implemented this scheme in a multiblock simulator. In light of the above modeling efforts, there is a need for a simulation method which more efficiently captures the effects of small scales on a coarse grid. Ideally, the method would be conservative and also treat tensor permeabilities correctly. Further, preferably the reconstructed finescale solution would satisfy the proper mass balance on the finescale. The present invention provides such a simulation method. SUMMARY OF THE INVENTION A multiscale finitevolume (MSFV) approach is taught for solving elliptic or parabolic problems such as those found in subsurface flow simulators. Advantages of the present MSFV method are that it fits nicely into a finitevolume framework, it allows for computing effective coarsescale transmissibilities, treats tensor permeabilities properly, and is conservative at both the coarse and finescales. The present method is computationally efficient relative to reservoir simulation now in use and is well suited for massive parallel computation. The present invention can be applied to 3D unstructured grids and also to multiphase flow. Further, the reconstructed finescale solution satisfies the proper mass balance on the finescale. A multiscale approach is described which results in effective transmissibilities for the coarsescale problem. Once the transmissibilities are constructed, the MSFV method uses a finitevolume scheme employing multipoint stencils for flux discretization. The approach is conservative and treats tensor permeabilities correctly. This method is easily applied using existing finitevolume codes, and once the transmissibilities are computed, the method is computationally very efficient. In computing the effective transmissibilities, closure assumptions are employed. A significant characteristic of the present multiscale method is that two sets of basis functions are employed. A first set of dual basis functions is computed to construct transmissibilities between coarse cells. A second set of locally computed finescale basis functions is utilized to reconstruct a finescale velocity field from a coarse scale solution. This second set of finescale basis functions is designed such that the reconstructed finescale velocity solution is fully consistent with the transmissibilities. Further, the solution satisfies the proper mass balance on the small scale. The MSFV method may be used in modeling a subsurface reservoir. A fine grid is first created defining a plurality of fine cells. A permeability field and other finescale properties are associated with the fine cells. Next, a coarse grid is created which defines a plurality of coarse cells having interfaces between the coarse cells. The coarse cells are ideally aggregates of the fine cells. A dual coarse grid is constructed defining a plurality of dual coarse control volumes. The dual coarse control volumes are ideally also aggregates of the fine cells. Boundaries surround the dual coarse control volumes. Dual basis functions are then calculated on the dual coarse control volumes by solving local elliptic or parabolic problems, preferably using boundary conditions obtained from solving reduced problems along the interfaces of the course cells. Fluxes, preferably integral fluxes, are then extracted across the interfaces of the coarse cells from the dual basis functions. These fluxes are assembled to obtain effective transmissibilities between coarse cells of the coarse cell grid. The transmissibilities can be used for coarse scale finite volume calculations. A finescale velocity field may be established. A finite volume method is used to calculate pressures in the coarse cells utilizing the transmissibilities between cells. Finescale basis functions are computed by solving local elliptic or parabolic flow problems on the coarse cells and by utilizing finescale fluxes across the interfaces of the coarse cells which are extracted from the dual basis functions. Finally, the finescale basis functions and the corresponding coarse cell pressures are combined to extract the finescale velocity field. A transport problem may be solved on the fine grid by using the finescale velocity field. Ideally, the transport problem is solved iteratively in two stages. In the first stage, a finescale velocity field is obtained from solving a pressure equation. In the second stage, the transport problem is solved on the fine cells using the finescale velocity field. A Schwartz overlap technique can be applied to solve the transport problem locally on each coarse cell with an implicit upwind scheme. A solution may be computed on the coarse cells at an incremental time and properties, such as a mobility coefficient, may be generated for the fine cells at the incremental time. If a predetermined condition is not met for all fine cells inside a dual coarse control volume, then the dual and finescale basis functions in that dual coarse control volume are reconstructed. BRIEF DESCRIPTION OF THE DRAWINGS These and other objects, features and advantages of the present invention will become better understood with regard to the following description, pending claims and accompanying drawings where: FIG. 1 illustrates a coarse 2D grid of coarse cells with an overlying dual coarse grid including a dual coarse control volume and an underlying fine grid of fine cells; FIG. 2 illustrates a coarse grid including nine adjacent coarse cells (bold solid lines) with a corresponding overlying dual coarse grid (bold dashed lines) including dual coarse control volumes and an underlying fine grid (thin dotted lines) of fine cells; FIG. 3 shows flux contribution q(2) Aand q(2)B due to the pressure in a particular coarse cell 2; FIG. 4 is a flowchart describing the overall steps used in a preferred embodiment of a reservoir simulation which employs a multiscale finitevolume (MSFV) method made in accordance with this invention; FIG. 5 is a flowchart further detailing steps used to determine transmissibilities T between coarse cells; FIG. 6 is a flow chart further describing steps used to construct a set of finescale basis functions and to extract a finescale velocity field; FIG. 7 is a flowchart depicting coupling between pressure and the saturation equations which utilize an implicit solution scheme and wherein Hand £ are operators used to update total velocity and saturation, respectively, during a single time step; FIG. 8 is an illustration of the use of an adaptive scheme to selectively update basis functions; FIG. 9 is an illustration of a permeability field associated with a SPE 10 problem; FIGS. 10AB are illustrations of permeability fields of a top layer and a bottom layer of cells from the SPE 10 problem; FIGS. 11AB are illustrations of saturation fields of top layers of cells created using the MSFV method and FIG. 11C is an illustration of a saturation field computed by a conventional finescale reservoir simulator; FIGS. 12AB are illustrations of saturation fields of bottom layers of cells created using the MSFV method and FIG. 12C is an illustration of a saturation field computed by a conventional finescale reservoir computer; FIGS. 13AB are graphs of oil cut and oil recovery; FIG. 14 is an illustration of a 3D test case having a grid of 10 x 22 x 17 grid cells and including injector and producer wells; and FIG. 15 is a graph of oil cut and oil recovery. BEST MODES FOR CARRYING OUT THE INVENTION I. FLOW PROBLEM A. One Phase Flow Fluid flow in a porous media can be described by the elliptic problem (Equation Removed) where p is the pressure, 4 is the mobility coefficient (permeability, K, divided by fluid viscosity, µ) and Ωis a volume or region of a subsurface which is to be simulated. A source term /represents wells, and in the compressible case, time derivatives. Permeability heterogeneity is a dominant factor in dictating the flow behavior in natural porous formations. The heterogeneity of permeability K is usually represented as a complex multiscale function of space. Moreover, permeability K tends to be a highly discontinuous full tensor. Resolving the spatial correlation structures and capturing the variability of permeability requires a highly detailed reservoir description. The velocity u of fluid flow is related to the pressure field through Darcy*s law: (Equation Removed) On the boundary of a volume, 9Ω , the flux q = u • v is specified, where v is the boundary unit normal vector pointing outward. Equations (1) and (2) describe incompressible flow in a porous media. These equations apply for both single and multiphase flows when appropriate interpretations of the mobility coefficient A and velocity u are made. This elliptic problem is a simple, yet representative, description of the type of systems that should be handled efficiently by a subsurface flow simulator. Moreover, the ability to handle this limiting case of incompressible flow ensures that compressible systems can be treated as a subset. 6. Two Phase Flow The flow of two incompressible phases in a heterogeneous domain may be mathematically described by the following: (Equation Removed) on a volume Q, where p is the pressure, So,w are the saturations (the subscripts o and w stand for oil and water, respectively) with 0 source terms which represent the wells. The system assumes that capillary pressure and gravity are negligible. Equivalently, system (3) can be written as: (Equation Removed)and the total mobility (Equation Removed) where ki = krj I>,. for j e [o, w}. Equation (4) is known as the "pressure equation" and equation (5) as the "hyperbolic transport equation." Again, equations (4) and (5) are a representative description of the type of systems that should be handled efficiently by a subsurface flow simulator. Such flow simulators, and techniques employed to simulate flow, are well known to those skilled in the art and are described in publications such as Petroleum Reservoir Simulation. K. Aziz and A. Settari, Stanford Bookstore Custom Publishing, 1999. II. MULTISCALE FINITEVOLUME (MSFV) METHOD A. MSFV Method for One Phase Flow 1. FiniteVolume Method A cell centered finitevolume method will now be briefly described. To solve the problem of equation (1), the overall domain or volume Ω is partitioned into smaller volumes (Ω) A finitevolume solution then satisfies (Equation Removed)for each control volume Ω where v is the unit normal vector of the volume boundary dΩ, pointing outward. The challenge is to find a good approximation for u v at 5Ω . In general, the flux is expressed as: (Equation Removed)Equation (9) is a linear combination of the pressure values, p , in the volumes {Ω} of the domain dΩ.The total number of volumes is n and Tk denotes transmissibility between volumes {Ω}. By definition, the fluxes of equation (9) are continuous across the interfaces of the volumes {Ω}and, as a result, the finitevolume method is conservative. 2. Construction of the Effective Transmissibilities The MSFV method results in multipoint stencils for coarsescale fluxes. For the following description, an orthogonal 2D grid 20 of grid cells 22 is used, as shown in FIG. 1. An underlying fine grid 24 of fine grid cells 26 contains the finescale permeability K information. To compute the transmissibilities T between coarse grid cells 22, a dual coarse grid 30 of dual coarse control volumes 32 is used. A control volume 32 of the dual grid 30, Ω, is constructed by connecting the midpoints of four adjacent coarse grid cells 22. To relate the fluxes across the coarse grid cell interfaces 34 which lie inside a particular control volume 32, or Ω, to the finitevolume pressures pk(k = l,4) in the four adjacent coarse grid cells 22, a local elliptical problem in the preferred embodiment is defined as (Equation Removed) For one skilled in the art, the method can easily be adapted to use a local parabolic problem. For an elliptic problem, Dirichlet or Neumann boundary conditions are to be specified on boundary dΩ Ideally, the imposed boundary conditions should approximate the true flow conditions experienced by the subdomain in the full system. These boundary conditions can be time and flow dependent. Since the subdomain is embedded in the whole system, Wallstrom et al. (T.C. Wallstrom, T.Y. Hou, M.A. Christie, L.J. Durlofsky, and D.H. Sharp, Application of a new twophase upscaling technique to realistic reservoir cross sections, SPE 51939, presented at the SPE Symposium on Reservoir Simulation, Houston, 1999) found that a constant pressure condition at the subdomain boundary tends to overestimate flow contributions from high permeability areas. If the correlation length of permeability is not much larger than the grid size, the flow contribution from high permeability areas is not proportional to the nominal permeability ratio. The transmissibility between two cells is a harmonic mean that is closer to the lower permeability. As a result, uniform flux conditions along the boundary often yield much better numerical results for a subdomain problem than linear or constant pressure conditions. Hou and Wu (T. Hou and W.H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comp. Phys, 134:169189, 1997) also proposed solving a reduced, problem to specify the boundary conditions for the local problem. The subscript t denotes the component parallel to the boundary of the dual coarse control volume 32 ora5 . For equation (1 1 ) and for the following part of this specification, Einstein summation convention will be used. The elliptic problem on a control volume D with boundary conditions of equation (11 ) on d& can be solved by any appropriate numerical method. In order to obtain a pressure solution that depends linearly on the pressures pk(j = 1,4), this preferred embodiment solves four elliptic problems, one for each cellcenter pressure. For instance, to get the solution for the pressure pl in the coarse grid cell having node 1 at its center, pk = Ω1k is set. The four solutions provide the dual basis functions O (k = 1,4) in control volume n , and the pressure solution of the local elliptic problem in a control volume Ω is the linear combination (Equation Removed Accordingly, the flux q across the grid cell interfaces can be written as a linear combination (Equation Removed) where qk(k 1,4) are the flux contributions from the corresponding dual basis functions, given all Ω (k = 1,4) from all control volumes Ω. The effective transmissibilities Tare computed, which can be used for finitevolume simulations, by assembling the flux contributions, in the preferred embodiment integral flux contributions across the cell interfaces 34. Note that the domain £2 can have any finescale distribution of mobility coefficients A. Of course the boundary condition given by equation (11) is an approximation that allows one to decouple the local problems. The MSFV and global finescale solutions are identical, only if equation (11) happens to capture the exact finescale pressure solution. However, numerical experiments have been performed which indicate that equation (11) is an excellent approximation of the boundary condition. Although the MSFV approach is a finitevolume method, it resembles the multiscale finiteelement method of Wu and Hou, briefly mentioned above. The construction of the dual basis functions is similar, though in the present MSFV method they are represented on the dual coarse grid rather than on the boundary of a finite element. A significant difference is that the present MSFV method is a cellcentered finitevolume method and is conservative. On the other hand, the mass matrix in the multiscale finiteelement method is constructed based on a variational principle and does not ensure local conservation. In the next section, the importance is illustrated of a finescale velocity field that is conservative. 3. Reconstruction of a Conservative FineScale Velocity Field Fluxes across the coarse cell interfaces 34 can be accurately computed by multiscale transmissibilities T. In some cases, it is interesting to accurately represent the smallscale velocities u (e.g., to predict the distribution of solute transported by a fluid). A straightforward approach might appear to be to use the dual basis functions Ω of equation (12). However, then the reconstructed finescale velocity field is, in general, discontinuous at the cell interfaces of the dual grid 30. Therefore, large errors can occur in the divergence field, and local mass balance is violated. Note that mass conservation is always satisfied for the coarse solution using the present MSFV method. The construction of a second set of local finescale basis functions Ω will now be described which is fully consistent with the fluxes q across the cell interfaces given by the dual basis functions Ω. This second set of finescale basis functions Ω allows a conservative finescale velocity field to be reconstructed. FIG. 2 shows a coarse grid 20 with nine adjacent grid cells 22 and a corresponding dual grid 30 of dual coarse control volumes 32 or Ω. For indexing purposes, these particular cells and corresponding dual volumes shall now be identified with numerals "19" and letters "AD" at their respective centers. Also shown is the underlying fine grid 24 of fine grid cells 26. The coarse grid, having the nine adjacent coarse cells 19, is shown in bold solid lines. The corresponding dual grid 30 of dual coarse control volumes AD are depicted with bold dashed lines. The underlying fine grid 24 of fine grid cells 26 is shown with thin dotted lines. To explain the reconstruction of the finescale velocity, the mass balance of the center grid cell 5 is examined. The coarse scale pressure solution, together with the dual basis functions Ω, provides the finescale fluxes q across the interfaces of coarse cell 5. To obtain a proper representation of the finescale velocity field in coarse cell 5, (i) the finescale fluxes across an interface of coarse cell 5 must match, and (ii) the divergence of the finescale velocity field within the coarse volume satisfies where QS is the coarse grid cell 5. The finescale flux q across the boundary of grid cell 5 depends on the coarse pressure solutions in grid cells 19. Therefore, the finescale velocity field within coarse grid cell 5 can be expressed as a superposition of finescale basis functions the help of FIG. 3, which depicts the needed dual coarse control volumes, the needed dual coarse control volumes, the construction the needed dual coarse control volumes, the construction of the finescale the needed dual coarse control volumes, the construction of the finescale the needed dual coarse control volumes, the construction of the finescaletion in construction of the finescale basis functions Φ' will be described. Each coarse cell pressure p(i 1,9) contributes to the finescale flux q . For example, let the contribution of the pressure in ceil 2 to the flux q in grid cell 5 beg(2). Note that q(2) is composed of contributions q(2)A and q(2)B coming from the dual basis functions associated with node 2 of volume A and volume B, respectively. To compute the finescale basis function Vz associated with the pressure in a coarse cell i, pJ = 6ij is set, and the pressure field is constructed according to the following equation. (Equation Removed) The finescale fluxes q are computed from the pressure field. These fluxes provide the proper boundary condition for computing the finescale basis function Φ. To solve the elliptic problem (Equation Removed) with the boundary conditions described above, solvability must be ensured. This is achieved by setting which is an equally distributed source term within Φ5. Finally, the solution of the elliptic problem, (16) and (17), is the finescale basis function & for coarse cell 5 associated with the pressure in volume /.The smallscale velocity field is extracted from the superposition (Equation Removed) For incompressible flow, this velocity field is divergence free everywhere. Computing the finescale basis functions Φrequires solving nine small elliptic problems, which are of the same size as those for the transmissibility calculations. Note that this step is a preprocessing step and has to be done only once. Furthermore, the construction of the finescale basis functions Φis independent and therefore well suited for parallel computation. The reconstruction of the finescale velocity field is a simple superposition and is ideally performed only in regions of interest. Alternatively, a conservative finescale velocity field may also be constructed directly in place. This construction may be performed as follows: (i) compute the finescale fluxes across the coarse cell interfaces using the dual basis functions with the pressures for the coarse cells; (ii) solve a pressure equation on each of the coarse cells using the finescale fluxes computed in step (i) as boundary conditions to obtain finescale pressures; (iii) compute the finescale velocity field from Darc/s law using the finescale pressures obtained in step (ii) with the underlying finescale permeability. The pressure solution of step (ii) may be performed on a system with larger support (e.g., by oversampling around the coarse cell). III. IMPLEMENTATION OF THE MSFV METHOD FIG. 4 is a flow chart summarizing the steps employed in a preferred embodiment in simulating a reservoir using the MSFV algorithm of this invention. The MSFV algorithm consists of six major steps: A. compute transmissibilities T for coarsescale fluxes (step 100); B. construct finescale basis functions (step 200); C. compute a coarse solution at a new time level; (step 300); D. reconstructs the finescale velocity field in regions of interest (step 400); E. solve transport equations (step 500); and F. recomputes transmissibilities and also the finescale basis functions in regions where the total mobility has changed more than a predetermined amount (step 600). Steps AD describes a twoscale approach. The methodology can be applied recursively with successive levels of coarsening. In cases of extremely fine resolution, this multilevel approach should yield scalable solutions. Parts E and F account for transport and mobility changes due to evolving phases and will be described in more detail below. A. Computing Transmissibilities for CoarseScale Fluxes  Step 100 The transmissibility calculations can be done in a stand alone module (Tmodule) and are well suited for parallel computation. The transmissibilities Tcan be written to a file for use by any finitevolume simulator that can handle multipoint flux discretization. Referring now to FIG. 5, a flowchart describes the individual steps which are undertaken to compute the transmissibilities T for a coarse scale model. First, a finescale grid having fine cells with an associated permeability field K are created (step 110). Next, a coarse grid, having coarse cells corresponding to the finescale grid, is created (step 120). The fine and coarse grids are then passed into a transmissibility or 7module. Dual coarse control volumes 5 are constructed (step 130), one for each node of the coarse grid. For each dual coarse control volume Q, dual or coarse scale basis functions O^ are constructed (step 140) by solving local elliptic problems (equation (10)) for each volume Q. This local elliptic problem, as described in section II.A.2 above, and the permeability field K associated with the fine grid are used and the boundary conditions corresponding to equation (11) are utilized (step 135) in solving the elliptic problem. In cases where the fine and coarse grids are nonconforming (e.g., if unstructured grids are used), oversampling may be applied. Finally, the integral coarse scale fluxes q across the interfaces of the coarse cells are extracted (step 150) from the dual basis functions Φ These integral coarse scale fluxes q are then assembled (step 160) to obtain MSFVtransmissibilities T between grid cells of the coarse grid. The computation of transmissibilities 7 can be viewed as an upscaling procedure. That is, the constructed coarse pressure solutions are designed to account for, in some manner, the finescale description of the permeability K in the original finescale grid model. Thus, part A step 100  computing transmissibilities, is preferably a separate preprocessing step used to coarsen the original finescale model to a size manageable by a conventional reservoir simulator. These transmissibilities T may be written to a file for later use. A finitevolume simulator that can handle multipoint flux discretization can then use these transmissibilities T. 6. Construction of FineScale Basis Function and Finescale Velocity Field Step 200 FIG. 6 is a flowchart describing the steps taken to construct a set of finescale basis functions Φ which can be isolated in a separate finescale basis function Φ module. These finescale basis functions Φ can then be used to create a finescale velocity field. This module is only necessary if there is an interest in reconstructing the finescale velocity field from the coarse pressure solution. As described in Section II.A.3 above, if the original dual basis functions Φ are used in reconstructing the finescale velocity field, large mass balance errors can occur. Here, steps are described to compute the finescale basis functions Φ, which can be used to reconstruct a conservative finescale velocity field. The procedure (step 200) of FIG. 4 follows the description of Section II.A.3 and has to be performed only once at the beginning of a simulation and is well suited for parallel computation. The finescale grid (step 210), with its corresponding permeability field K, the coarse grid (step 220), and the dual basis functions Φ (step 230) are passed into a finescale basis function Φ . A pressure field is constructed from the coarse scale pressure solution and dual basis functions (step 250). The finescale fluxes for the coarse cells are then computed (step 260). For each control volume, elliptic problems are solved, using the finescale fluxes as boundary conditions, to determine finescale basis functions (step 270). The finescale velocity field can then be computed from the superposition of cell pressures and finescale basis functions. The results may then be output from the module. Alternatively, the finescale velocity field can be computed directly in place as has been described above in section II.A.3. In many cases, the finescale velocity field has to be reconstructed in certain regions only, as will be described in fuller detail below. Therefore, in order to save memory and computing time, one can think of an in situ computation of the finescale basis functions 3) , which, once computed, can be reused. C. Computation of the Coarse Solution at the New Time  Step 300 Step 300 can be performed by virtually any multipoint stencil finitevolume code by using the MSFVtransmissibilities 7for the flux calculation. These coarse fluxes effectively capture the largescale behavior of the solution without resolving the small scales. D. Reconstruction of the FineScale Velocity Field  Step 400 Step 400 is straightforward. Reconstruction of the finescale velocity field in regions of interest is achieved by superposition of the finescale basis FIG. 6. Alternatively, the finescale velocity field can be computed directly in functions Φas described in section II.A.3, step B above and as shown in place as described above in section II.A.3. Of course, many variations of the MSFV method can be devised. It may be advantageous; however, that construction of the transmissibilities T and finescale basis functions Φ can be done in modules separate from the simulator. E. Solving Pressure and Transport Equations 1. Numerical solution algorithm  explicit solution Multiphase flow problems may be solved in two stages. First, the total velocity field is obtained from solving the pressure equation (4), and then the hyperbolic transport equation (5) is solved. To solve the pressure equation, the MSFV method, which has been described above is used. The difference from single phase flow is that in this case the mobility term reflects the total mobility of both phases, and then the obtained velocity field u is the total velocity in the domain. The reconstructed finescale velocity field u is then used to solve the transport equation on the fine grid. The values of k0iW are taken from the upwind direction; time integration may be obtained using a backward Euler scheme. Note that, in general, the dual and finescale basis functions (6.O) must be recomputed each time step due to changes in the saturation (mobility) field. 2. Numerical Solution Algorithm  Implicit Coupling In the preferred embodiment of this invention, the MSFV method utilizes an algorithm with implicit calculations. The multiphase flow problem is solved iteratively in two stages. See FIG. 7 for a diagram of this method illustrating the coupling between the pressure and saturation equations. First, in each Newton step, a saturation field S is established  either initial input or through iteration (step 51 0). Next, a pressure equation (see equation (19) below) is solved (step 520) using the MSFV techniques described above to obtain (step 530) the total velocity field. Then a transport equation (see equation (20) below) is solved (step 540) on the fine grid by using the reconstructed finescale velocity field u. In this solution, a Schwarz overlap technique is applied, i.e., the transport problem is solved locally on each coarse volume with an implicit upwind scheme, where the saturation values from the neighboring coarse volumes at the previous iteration level are used for the boundary conditions. Once the Schwarz overlap scheme has converged (steps 550, 560) — for hyperbolic systems this method is very efficient — the new saturation distribution determines the new total mobility field for the pressure problem of the next Newton iteration. Note that, in general, some of the basis functions have to be recomputed each iteration. The superscripts n and v denote the old time and iteration levels, respectively. Saturation is represented byS , the total velocity field by u, the computation of the velocity by the operator!! , and the computation of the saturation by . The new pressure field pv+1 is obtained by solving (Equation Removed) from which the new velocity field u v+1 is computed. The new saturation field is obtained by solving (Equation Removed) F. Recomputing Transmissibilities and FineScale Basis Functions  Adaptive Scheme The most expensive part of the MSFV algorithm for multiphase flow is the reconstruction of the coarse scale and finescale basis functions (Φ ). Therefore, to obtain higher efficiency, it is desirable to recompute the basis functions only where it is absolutely necessary. An adaptive scheme can be used to update these basis functions. In the preferred exemplary embodiment, if the condition (Equation Removed) is not fulfilled (the superscripts n and «l denote the previous two time steps and s^ is a defined value) for all fine cells inside a coarse dual volume, then the dual basis functions of that control volume have to be reconstructed. Note that condition (23) is true if Φ changes by a factor which is larger than £%) and smaller than l +Φ. An illustration of this scheme is shown in FIG. 8, where the fine and the coarse grid cells are drawn with thin and bold lines, respectively. The black squares represent the fine cells in which condition (23) is not fulfilled. The squares with bold dashed lines are the control volumes for which the dual basis functions have to be reconstructed. The shaded regions represent the coarse cells for which the finescale basis functions have to be updated. In the schematic 2D example of FIG. 8, only 20 of 196 total dual basis functions and 117 of 324 total finescale basis functions have to be reconstructed. Of course, these numbers depend heavily on the defined thresholΦ. In general, a smaller threshold triggers more fine volumes, and as a consequence more basis functions are recomputed each time step. For a wide variety of test cases, it has been found that taking e% to be IV. NUMERICAL RESULTS This MSFV method, combined the implicit coupling scheme shown in FIG. 7, has been tested for two phase flow (µ0//zwsio) in a stiff 3D model with more than 140,000 fine cells. It has been demonstrated that the multiscale results are in excellent agreement with the finescale solution. Moreover, the MSFV method has proven to be approximately 27 times more efficient than the established oil reservoir simulator Chears. However, in many cases the computational efficiency is compromised due to the time step size restrictions inherent for IMPES schemes. This problem may be resolved by applying the fully implicit MSFV method, which was described in the previous section. Here numerical studies show the following: (1) The results obtained with the implicit MSFV method are in excellent agreement with the finescale results. (2) The results obtained with the implicit MSFV method are not very sensitive to the choice of the coarse grid. (3) The implicit MSFV for two phase flow overcomes the time step size restriction and therefore very large time steps can be applied. (4) The results obtained with the implicit MSFV method are, to a large extent, insensitive to the time step size. (5) The implicit MSFV method is very efficient. For the finescale comparison runs, the established reservoir simulator Chears was used. The efficiency of both the implicit MSFV method and the finescale reservoir simulator depends on the choice of various parameter settings which were not fully optimized. A. Test Case To study the accuracy and efficiency of the fully implicit MSFV algorithm, 2D and 3D test cases with uniformly spaced orthogonal 60 x 220 and 60 x 220 x 85 grids were used. The 3D grid and permeability field are the same as for the SPE 10 test case, which is regarded as being extremely difficult for reservoir simulators. While this 3D test case is used for computational efficiency assessment, the 2D test cases, which consist of top and bottom layers, serves to study the accuracy of the MSFV method. FIG. 9 illustrates the 3D test case posed by the permeability field of the SPE 10 problem. The darker areas indicate lower permeability. An injector well is placed in the center of the field and four producers in the corners. These well locations are used for all of the following studies. The reservoir is initially filled with oil and /// B. 2D Simulation of the Top and Bottom Layers The MSFV simulator used lacked a sophisticated well model. That is, wells are modeled by defining the total rates for each perforated coarse volume. Therefore, in order to make accuracy comparisons between MSFV and finescale (Chears reservoir simulator) results, each finescale volume inside each perforated coarse volume becomes a well in the Chears runs. For large 3D models, this poses a technical problem since Chears reservoir simulator is not designed to handle an arbitrary large number of individual wells. For this reason, it was determined to do an accuracy assessment in 2D, i.e., with the top and the bottom layers of the 3D model. These two layers, for which the permeability fields are shown in FIGS. 10A and 10B, are representative for the two characteristically different regions of the full model. MSFV simulations were performed with uniformly spaced 10 x 22 and 20 x 44 coarse grids. The results were compared with the finescale solution on a 60 x 220 grid. As in the full 3D test case, there are four producers at the corners which are distributed over an area of 6 x 10 finescale volumes. The injector is located in the center of the domain and is distributed over an area of 12 x 12 finescale volumes. The rates are the same for all finescale volumes (positive for the producer volumes and negative for the injector volumes). FIGS. 11AC and 12AC show the permeability fields of the respective top and the bottom layers. The black is indicative of low permeability. These two layers are representative for the two characteristically different regions of the full 3D model. FIGS. 11AC and 12AC show the computed saturation fields after 0.0933 PVI (pore volume injected) for the top and the. bottom layers, respectively. While FIGS. 11C and 12C show the finescale reference solutions, FIGS 11A and 11B and 12A and 12B show the MSFV results for 10 x 22 and 20 x 44 coarse grids, respectively. For both layers, it can be observed that the agreement is excellent and that the multiscale method is hardly sensitive to the choice of the coarse grid. A more quantitative comparison is shown in FIGS. 13A and 13B where the finescale and multiscale oil cut and oil recovery curves are plotted. Considering the difficulty of these test problems and the fact that two independently implemented simulators are used for the comparisons, this agreement is quite good. In the following studies, it will be demonstrated that for a model with 1,122,000 cells, the MSFV method is significantly more efficient than finescale simulations and the results remain accurate for very large time steps. C. 3D Simulations While 2D studies are appropriate to study the accuracy of the implicit MSFV method, large and stiff 3D computations are required for a meaningful efficiency assessment. A 3D test case was employed as described above. A coarse 10 x 22 x 17 grid, shown in FIG. 14, was used and 0.5 pore volumes were injected. Opposed to the MSFV runs, the wells for the CHEARS simulations were defined on the finescale. Table 1 below shows CPU time and required number of times steps for the CHEARS simulation and two MSFV runs. TABLE 1 Efficiency Comparison Between Msfv And FineScale Simulations (Table Removed) While Chears uses a control algorithm, the time step size in the multiscale simulations was fixed. It is due to the size and stiffness of the problem that much smaller time steps have to be applied for a successful Chears simulation. The table shows that the implicit MSFV method can compute the solution approximately 27 times faster than CHEARS. FIG. 15 shows the oil cut and recovery curves obtained with multiscale simulations using 50 and 200 time steps. The close agreement between the results confirms that the method is very robust in respect to time step size. Since the cost for MSFV simulation scales almost linearly with the problem size and since the dual and finescale basis function can be computed independently, the method is ideally suited for massive parallel computations and huge problems. While in the foregoing specification this invention has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, it will be apparent to those skilled in the art that the invention is susceptible to alteration and that certain other details described herein can vary considerably without departing from the basic principles of the invention. WE CLAIM: 1. A multiscale finitevolume method for use in modeling a subsurface reservoir comprising: (a) creating a fine grid defining a plurality of fine cells and having a permeability field associated with the fine cells; (b) creating a coarse grid defining a plurality of coarse cells having interfaces between the coarse cells, the coarse cells being aggregates of the fine cells; (c) creating a dual coarse grid defining a plurality of dual coarse control volumes, the dual coarse control volumes being aggregates of the fine cells and having boundaries bounding the dual coarse control volumes; (d) calculating dual basis functions on the dual coarse control volumes by solving local elliptic or parabolic problems; (e) extracting fluxes across the interfaces of the coarse cells from the dual basis functions; (f) assembling the fluxes to calculate effective transmissibilities between coarse cells; (g) calculating pressure in the coarse cells using a finite volume method and utilizing the effective transmissibilities between coarse cells; and (h) computing a finescale velocity field. 2. The method of claim 1 wherein: the finescale velocity field is computed directly in place. 3. The method of claim 2 wherein: the step of computing a finescale velocity field directly in place includes: (i) computing finescale fluxes across the coarse cell interfaces using the dual basis functions with the pressures the coarse cells; (ii) solving a pressure equation on each of the coarse cells using the finescale fluxes computed in step (i) as boundary conditions to obtain finescale pressures; and (iii) computing the finescale velocity field from Darcy's law using the finescale pressure obtained in step (ii). 4. The method of claim 3 wherein: the solving of the pressure equation of step (ii) to obtain finescale pressures is performed on a system with larger support. 5. The method of claim 4 wherein: the solving of the pressure equation of step (ii) may be performed by oversampling around the coarse cells. 6. The method of claim 3 wherein: the solving of the pressure equation of step (ii) may be performed by oversampling around the coarse cells. 

4588delnp2007Abstract(11022013).pdf
4588delnp2007Assignment(19032012).pdf
4588delnp2007Claims (11022013).pdf
4588delnp2007Correspondence Others(11022013).pdf
4588delnp2007Correspondence Others(13082012).pdf
4588delnp2007Correspondence Others(19032012).pdf
4588delnp2007correspondenceothers.pdf
4588delnp2007description (complete).pdf
4588delnp2007Drawings(11022013).pdf
4588delnp2007Form1(19032012).pdf
4588delnp2007Form2(19032012).pdf
4588delnp2007Form3(13082012).pdf
4588delnp2007Form5(19032012).pdf
4588delnp2007GPA(19032012).pdf
4588delnp2007Petition137(13082012).pdf
Patent Number  258810  

Indian Patent Application Number  4588/DELNP/2007  
PG Journal Number  07/2014  
Publication Date  14Feb2014  
Grant Date  08Feb2014  
Date of Filing  15Jun2007  
Name of Patentee  CHEVRON U.S.A. INC.  
Applicant Address  6001 BOLLINGER CANYON ROAD, 3RD FLOOR, SAN RAMON, CA 94583, USA  
Inventors:


PCT International Classification Number  C06G 7/48  
PCT International Application Number  PCT/US2005/042632  
PCT International Filing date  20051121  
PCT Conventions:
