MX2007007052A - Finite volume method for coupled stress/fluid flow in a reservoir simulator - Google Patents

Abstract

A method for conducting a stress calculation is disclosed adapted for modeling a set of stresses and displacements in a reservoir, the method comprising:(a) building a reservoir model over a region of interest by gridding the region of interest, the grid being comprised of one or more cells and having nodes, each cell having a cell center;(b) interpolating unknown rock displacements in the region of interest from cell centers to grid nodes;(c) integrating over each cell to form a discrete system of equations;and (d) using the discrete system of equations to model the stresses and displacements in the reservoir.

Description

SYSTEM AND METHOD WITH FINITE VOLUME AND PROGRAM STORAGE DEVICE FOR LINEAR ELASTICITY INVOLVING TENSION AND FLOW COUPLED IN A SIMULATOR OF DEPOSIT BACKGROUND The subject matter of this specification refers to a new finite volume method (and associated program storage system and device) for a set of Linear Elasticity equations described below; and in particular to a finite volume approach for discretization of a set of Linear Elasticity equations in a general unstructured grid, in three dimensions. The flow equations in commercial deposit simulators are generally discretized using the volume or finite difference method, while for stress equations, it is more common to use the finite element approach. In a deposit simulator designed to include geomechanical effects, these two different types of equations must be solved with some degree of coupling. It is therefore natural to ask: if suitable finite volume methods can be derived for effort equations, so that fluid and effort flow models can share a common derivation and what are the relative merits of the finite element and finite volume methods for these coupled systems. In this specification, a finite volume discretization of a set of stress equations, as implemented in a commercial warehouse simulator, is presented as an option for coupling the effort with fluid flow. The method, as presented, is locally conservative and retains second-order accuracy in general three-dimensional grids. The imposition of various types of boundary conditions is discussed and the implementation of special features is described, such as faults, narrowing and local grid refinement. A comparison with other approaches is presented based on finite differences (see references 1 and 2 below) and finite elements (see reference 3 below). The relative accuracy, efficiency and robustness of these three different approaches are also discussed. COMPENDIUM One aspect of the "Finite Volume Method for Linear Elasticity" as described in this specification, involves a method to conduct an effort calculation adapted to model a set of efforts and displacements in a deposit, which it comprises: (a) constructing a deposit model on a region of interest by gridding the region of interest, the grid comprising one or more cells and having nodes, each cell having a cell center; (b) interpolate unknown rock displacements in the region of interest from cell centers to grid nodes; (c) integrate over each cell to form a discrete system of equations; and (d) use the discrete system of equations to model displacements and stresses in the deposit. Another aspect of the "Finite Volume Method for Linear Elasticity" as described in this specification, involves a machine-readable storage device that tangibly incorporates a set of instructions that are executed by machine to perform method steps to To conduct an effort calculation adapted to model a set of efforts and displacements in a deposit, the method steps comprise: (a) constructing a deposit model on a region of interest by gridding the region of interest, the grid comprises a or more cells and have nodes, each cell has a cell center; (b) interpolate unknown rock displacements in the region of interest from cell centers to grid nodes; (c) integrate over each cell to form a discrete system of equations; and (d) use the discrete system of equations to model displacements and stresses in the deposit. Another aspect of the "Finite Volume Method for Linear Elasticity" as described in this specification, involves a method to conduct an effort calculation to model a set of efforts and displacements in the tank, which comprises: first apparatus adapted to build a model of deposits in a region of interest by grid-forming the region of interest, the grid comprises one or more cells, and has nodes, each cell has a cell center; a second apparatus adapted to interpolate unknown rock shifts in the region of interest from cell centers to grid nodes; a third apparatus adapted to integrate on each cell to form a system of discrete equations; and a fourth apparatus adapted to model the deposit while the discrete system of equations is used to model displacements and stresses in the deposit. Another aspect of the "Finite Volume Method for Linear Elasticity" as described in this specification, involves a computer program adapted to be executed by a processor, the computer program when executed by a processor, practice a process adapted to conduct a calculation of effort to model a set of efforts and displacements in a deposit, the process comprises: (a) building a deposit model on a region of interest by grid-forming the region of interest, the grid it comprises one or more cells and has nodes, each cell has a cell center; (b) interpolating unknown rock shifts in the region of interest from cell centers to grid nodes; (c) integrate over each cell to form a discrete system of equations; and (d) use the discrete system of equations to model the efforts and displacements in the deposit. Additional scope of applicability will be apparent from the detailed description presented below. It should be understood, however, that the detailed description and specific examples set forth below are given by way of illustration only, since various changes and modifications within the spirit and scope of the "Finite Volume Method for Linear Elasticity" as described and claimed in this specification, will be apparent to a person with technical skill, from a reading of the following detailed description. BRIEF DESCRIPTION OF THE DRAWINGS A complete understanding will be obtained from the detailed description presented here below and the accompanying drawings which are given by way of illustration only and are not intended to be limiting in any proportion, and wherein: Figure 1 shows a seismic operation for producing a reduced seismic data output register, the seismic operation of Figure 1 includes a data reduction operation; Figure 2 illustrates a drilling operation to produce a log output log against depth or time or both of the well; Figure 3 illustrates a computer system for performing the data reduction operation of Figure 1; Figures 4 and 5 illustrate a workstation adapted to store a program or software "Flogrid" and a simulator software "Eclipse"; Figures 6 and 7 illustrate a more detailed construction of the "Flogrid" program of Figure 5 which is adapted to generate output data, for use by the "Eclipse" simulator program, the Eclipse simulator program includes a "Finite Volume Method for Linear Elasticity Equations "that are described in this specification; Figure 8 illustrates an example of a typical output display generated by the "Eclipse" simulator program of Figure 6, which is displayed on the 3D viewer of Figure 6; Figure 9 illustrates a prior art approach or method for performing deposit simulation that has been practiced by the deposit simulators of the prior art; Figure 10 illustrates the Eclipse simulator program of Figures 5 and 6, which includes the "Finite Volume Method for Linear Elasticity Equations"; as described in this specification; Figure 11 and illustrates a section of a Two-dimensional Finite Volume grid; Figures 12 and 13 illustrate grid and narrowing sections and Local Grid Refinement; Figure 14 illustrates integration on a control volume face; Figure 15 illustrates a contour in sliced plane through grating with local refinement; Figure 16 illustrates grid transformed with narrowing; Figure 17 illustrates a slice plane through grid with random distortion; Figure 18 illustrates a transformed grid with narrows, faults and local grid refinements; Figure 19 illustrates the rate of convergence in a sequence of grids of various types for a general known real solution; Figure 20 illustrates a table showing multiple grid interactions (upe time) for a sequence of meshes and transformations; Figure 21 illustrates the term of Effective Effort s ?? for the finite volume solution within a commercial simulator; Figure 22 illustrates the term of Effective Effort s ?? for the finite volume solution within the commercial simulator; Figure 23 illustrates the term of Effort Cash s ?? for finite volume solution within the commercial simulator with Local Grid Refinement around wells; DESCRIPTION Oil and gas are produced from underground rock formations. These rocks are porous, such as sponge, and are filled with fluid, usually water. This porous feature of the rocks is known as porosity. These rocks, besides being porous, have the ability to allow fluids to flow through pores, a characteristic measured by a property called permeability. When oil (or gas) is trapped in these formations, it may be possible to extract it when drilling wells that derive in the formation. As long as the pressure in the well is less than in the rock formation, the fluids contained in the pores will flow into the well. These fluids can flow naturally up the well to the surface or the upflow into the well may have to be assisted by pumps. The relative amounts of oil, gas and water that are produced on the surface will depend on the fraction of the rock pore space occupied by each type of fluid. Water is always present in the pores, but it will not flow unless its volume fraction exceeds a threshold value that varies from one type of rock to another. Similarly, oil and gas will only flow as long as their volume fractions exceed their own thresholds. The characteristics of the rock (including porosity and permeability) in an oil deposit vary enormously from one site to another. As a result, the relative amounts of oil, gas and water that can be produced will also vary from deposit to deposit. These variations make it difficult to simply forecast the amount of fluid and gases that a deposit will produce and the amount of resources it requires to produce from a particular deposit. However, interested parties to produce from a deposit require to project the production of the deposit with some precision in order to determine the feasibility of producing from that deposit. Therefore, in order to accurately predict production speeds of all the wells in a deposit, it is necessary to construct a detailed mathematical model of the geology and geometry of the deposits. A large amount of research has focused on the development of tools for deposit simulation. These tools include mathematical and computer models, which they describe and which are used to forecast the flow of multiple phases of oil and gas within a three-dimensional underground formation (a "field"). Deposit tools use data acquired empirically to describe a field. These data are combined with and manipulated by mathematical models whose output describes specified characteristics of the field as a future time and that in terms of quantities to be measured such as the production or injection speeds of individual wells and groups of wells, head pressure in the pipeline or at the bottom of the well in each well and the distribution of pressure and fluid phases within the reservoir. The mathematical model of a deposit typically It is done by dividing the deposit volume into a large number of interconnected cells and estimating the average permeability, porosity and other rock properties for each cell. This process uses seismic data, well logs and rock cores recovered when wells are drilled. The production of the deposit can then be mathematically modeled by numerically solving a system of three or more non-linear partial differential equations that describe the fluid flow in the reservoir. Computer analysis of production of an oil deposit, usually divided into two phases, coupling or historical correspondence and prediction. In the historical coupling phase, the past production behavior of the deposit and its wells is modeled repeatedly, starting with initial production and continuing until the current time. The first computer run, is based on a geological model as described above. After each run, the computer results are compared in detail with data collected in the oil field during the entire production period. Geo-scientists modify the geological model in deposit, based on the differences between current and calculated production performance and re-run the computer model. East process continues, until the mathematical deposit model behaves like the real oil deposit. Once a convenient historical correspondence has been obtained, the production of the oil deposit can be forecasted far into the future (sometimes for up to 50 years). Oil recovery can be maximized and production costs minimized by comparing many alternate operating plans, each requiring a new run of the computer model. After a field development plan is put into action, the deposit model can be periodically re-executed and also adjusted to improve its ability to match newly collected production data. When sufficient data is obtained regarding the deposit, the characteristics of a deposit can be mathematically moderated to predict production speeds from wells in that deposit. The gross characteristics of the field include the porosity and permeability of the reservoir rocks, the thickness of the geological zones, the location and characteristics of the geological faults, relative permeability and capillary pressure functions and the characteristics of reservoir fluids such as density, viscosity and relationships phase balance. From these data, a set of continuous partial differential equations (PDEs = partial differential equations) is generated, which describes the behavior of the field as a function of production and time parameters. These production parameters include the location of wells, the characteristics of well completions and the operational restrictions applied to the wells. Operational constraints may include such as the production rate of a particular fluid phase, the bottom hole pressure, the pipe head pressure or the combined flow expense of a group of wells. These restrictions can be applied directly by data or by another simulator that models the flow of fluids in the surface equipment used to transport the fluids produced from or injected into the wells. However, because only the simplest system of PDEs can be solved using classical or closed-form techniques (for example a homogeneous field that has circular borders), the PDEs of a model become a set of non-linear approximations that then they are solved numerically. One approximate technique is the finite difference method. In the finite difference method, deposit PDEs are converted into a series of different ratios that they divide a deposit into a collection of discrete three-dimensional cells, which are then resolved by discrete times to determine (or forecast) the value of deposit characteristics such as pressure, permeability, fluid fractions and at a later time. As oil and gas are extracted from a reservoir, through the wells, the pressure in reservoirs is reduced, and there is a tendency for the porous rock to be compacted in the formation. For strong rocks, this effect may not have a significant effect on the simulation and does not need to be taken into consideration when solving the reservoir flow equations. For weaker rocks, the effect is more significant and ideally we want to model the movement of the reservoir rock by solving the linear elasticity equations that describe its displacements and the forces present in the rock. These equations are coupled with the reservoir flow equations for a pressure term in the stress tensor. There is also an implicit coupling by virtue of the fact that the grid cells used to solve the flow equations should ideally move according to the displacements obtained from solving the elasticity equations. In practice, for small relative displacements, this effect can be reasonably approximated by a modification to the porosity. Inside the computerized "deposit simulator", the performance of deposits is modeled in discrete increments of time. Each time stage so named advances the solution from a previous point in time, where all the variables are known, to a future point in time where all the variables are unknown. This process is repeated until the entire time period of interest has been modeled. Within each time step, it is necessary to solve a large system of nonlinear equations that model the fluid flow, from cell to cell and through the wells. (With the current technology, it is possible to include several million cells in the deposit model). Solutions to the system of nonlinear equations are obtained by Newton's iteration. In each iteration the system of nonlinear equations is approximated by a system of linear equations, which must be solved by yet another interactive procedure. A similar "deposit simulator" is the "eclipse" deposit simulator that is owned and operated by Schlumberger Technology Corporation of Houston, Tex. The "eclipse" simulator program receives data from output of the grid training program for simulation of "flogrid" and in response to this, the simulator program "eclipse" generates a set of simulation results that are displayed in a 3D viewer. The "flogrid" simulation grid formation program is described in US Pat. No. 6,106,561 issued to Farmer, the description of which is hereby incorporated by reference in this specification. As illustrated in Figure 10, the "eclipse" simulator program includes an "infinite volume method for linear elasticity equations" program (including its associated program and storage device). The "infinite volume method for linear elasticity equations" of Figure 10 and associated program storage system and device) describes in detail a possible implementation of an infinite volume method for linear elasticity, indicating its convenience in non-threaded meshes. Structured and knitted corner and shows your order of precision. This specification includes: (1) a background layout with reference to Figures 1-9 that provides background information regarding the performance of a seismic operation and a well registration operation adapted to generate registration data of well and seismic, well and seismic log data are provided as power data to a workstation that stores a "flogrid" simulation grid formation program and an "eclipse" simulator program, and (2) a description of the "eclipse" simulator program also includes the "Finite Volume Method" program "for linear elasticity equations" (and associated program storage system and device) discussed below with reference to Figures 10-23 which represents an implementation possible from an "infinite volume method for linear elasticity". With reference to Figure 1, a method and apparatus for performing a seismic operation are illustrated. During a seismic operation, a source of acoustic energy or sound vibrations 10, such as an explosive energy source 10, produces a plurality of sound vibrations. In Figure 1, a similar sound vibration 12 is reflected from a plurality of horizons 14 in a terrestrial formation 16. The sound vibrations or sound 12 are received in a plurality of geophonic receivers 18 located on the surface of the earth and the geophones 18 produce electrical output signals, referred to as "received data" 20 in Figure 1, in response to the sound or vibrations received 12 representative of different parameters (such as amplitude and / or frequency) of the sound vibration (s) 12. The "received data" 20 are provided as "feeding data" to a computer 22a of a registration truck 22 and in In response to the "feed data", the log truck computer 22a generates a "seismic data output log" 24. Subsequently in the processing of the seismic data output log 24, these seismic data are subjected to "reduction of data "30 on a large computer and a" reduced seismic data output register "24a, is generated from this data reduction operation 30. Referring to Figure 2, a well registration operation is illustrated. During the well registration operation, a well registration tool 34 is lowered into the land formation 16 of Figure 1, which is penetrated by a bore 36. In response to the well registration operation, registration data of well 38 are generated from the well registration tool 34, the well registration data 38 is provided as "power data" to a computer 40a of a well log 40 truck. In response to the log data of well 38, the 40a well log truck computer produces a "Well log output log" 42. With reference to Figure 3, the registration register of seismic data 24 of Figure 1 is provided as "power data" to a large computer 30, where the data reduction operation 30 of Figure 1 is performed. A large computer processor 30a will execute a data reduction program 30b stored in a large 30b computer storage. When the execution of the data reduction program 30b is completed, the reduced seismic data output register 24a of FIGS. 1 and 3 is generated. With reference to FIGS. 4 and 5, a workstation 44 is illustrated in FIG. 4. A storage medium 46 such as a CD-ROM 46 stores the program, and that program can be loaded on the workstation 44, for storage in the workstation memory. In Figure 5, the workstation 44 includes a workstation memory 44a, the program stored in the storage medium (CD-ROM) 46 is loaded into the workstation 44 and stored in the memory of the workstation. work 44a. A workstation processor 44d will execute the program stored in the workstation memory 44a, in response to certain power data that is provide the workstation processor 44d, and then the processor 44d will display or record the results of that processing in the "3D viewer or display or recorder" of the workstation 44e. The power data, which is provided to the workstation 44 in Figure 5, includes the log output from the well log 42 and the reduced seismic data log output 24a. The "log out of the well log" 42 represents the well log data generated during the logging operation in a landform of Figure 2, and the "log of reduced seismic data" 24a represents seismic data reduced in data generated by the large computer 30 in Figure 3, in response to the seismic operation illustrated in Figure 1. In Figure 5, the program stored in the storage medium (CD-ROM) 46 in Figure 5 includes a "flogrid" program 46a and an "eclipse" simulator program 46b. When the storage medium (CD-ROM) 46 is inserted into the workstation 44 in Figure 5, the "flogrid" program 46a and the "eclipse" simulator program 46b stored on CD-ROM 46, both are loaded into the workstation 44 and stored in the workstation memory 44a. The program "flogrid" 46a and the simulator program "eclipse" 46b are owned and operated by Schlumberger Technology Corporation of Houston, Tex. Program "flogrid" 46a is described in the patent of the US. No. 6,106,561 granted to Farmer under the title "Simulation Gridding Method and Apparatus including a Structured A real Gridder Adapted for use by a Reservoir Simulator", the description of which is incorporated by reference in this specification. When the workstation processor 44d executes the flogid program 46a and the eclipse simulator program 46b, the "eclipse" simulator program 46b responds to a more accurate set of grid information information associated with a respective set of grid blocks of a structured simulation grid generated by the "flogrid" program 46a, by additionally generating a set of more precise simulation results that are respectively associated with the set of grid blocks of the simulation grid. These simulation results are displayed on the 3D viewer 44e of Figure 5 and can be recorded on a recorder 44e. With reference to Figures 6 and 7, initially with reference to Figure 6, the flogrid program 46a and the eclipse simulator program 46b are illustrated as stored in the workstation memory 44a of Figure 5. In addition, the "results of simulation "48 that are sent out of the program eclipse 46b simulator, are illustrated received by and displayed on the 3d 44e viewer. The flogrid program 46a includes a deposit data store, a deposit frame, a structured grid generator, an unstructured grid generator, and a scale amplifier, all of which are fully discussed in U.S. Pat. No. 6,106,561 granted to Farmer previously referred to, the description of which has already been incorporated by reference in this specification. In Figure 6, a set of "grids and simulation properties associated with the grids" 47, generated by the scale amplifier and the unstructured grid generator "petragrid" are received by the eclipse simulator program 46b. In response, the eclipse simulator program 46b generates a "set of simulation results associated respectively by a set of grid blocks of the simulated grids" 48, and the simulation results and the associated grid blocks 48 are displayed on the 3D viewer 44e. The unstructured grid former "Petragrid" is described in U.S. Pat. Nos. 6,018,497 and 6,078,869, the descriptions of which are incorporated by reference in this specification. In Figure 7, the flogrid 46a program generates an output data set 47, comprising a plurality of grid cells and certain properties associated with grid structures. This output data 47 is provided as power data to the eclipse simulator program 46b. Some other programs 49 provide other power data to the eclipse 46b simulator program. In response to the output data 47 (comprising a terrestrial array consisting of a grid, including a plurality of grid cells and certain properties associated with each of the grid cells), as well as the other output data of the other programs 49, the eclipse simulator program 46b generates a set of "simulation results" 48, the simulation results 48 include the plurality of grid cells and a plurality of simulation results associated, respectively, with the plurality of grid cells. The aforementioned plurality of grid cells and the plurality of simulation results associated, respectively with the plurality of grid cells, are displayed in the 3D viewer 44e of Figures 6 and 7. With reference to Figure 8, an example of the simulation results 48 (ie the "plurality of grid cells" and the plurality of results of associated simulations, respectively with the plurality of grid cells ") that are displayed in the 3D viewer 44e of Figures 5, 6 and 7, is illustrated in Figure 8. The following paragraphs will pre the eclipse simulator program 46b, of the Figures 5, 6 and 7, wherein the eclipse simulator program 46b further includes a detailed description of a "finite volume method" (and associated program and storage device) "for linear elasticity equations", as illustrated in Figure 10. In particular, the following paragraphs will pre a detailed description of a finite-volume approach to discretization of the linear elasticity equations in a general unstructured grid, in three dimensions, as shown in Figure 11. With reference to Figure 10, a general operational profile of a reservoir simulator of the prior art, is discussed below with reference to Figure 9. In Figure 9, the data from reservoir 42 and 24a of Figure 5 and the rock core data is used to describe a computational grid and the properties of the reservoir rocks. These data are combined with data referring to the physical properties of the fluids contained in the deposit, the data combined are used to calculate the initial distributions of pressure and fluid saturations (volume fractions) as well as the composition of each fluid phase, block 50 in Figure 9. Data of variation over time, such as locations and characteristics of wells, production and flow controls of injection flow, and simulator control information, are read from a database, block 52. Using the current pressure, saturation and fluid compositions for each grid cell, the partial differential equations that describe Mass balances are approximated by finite differences in block 54, which results in two or more non-linear algebraic equations for each grid cell. Also, in block 54, these non-linear equations are linearized by Newton's method. In block 56, the resulting system of linear equations is solved iteratively, using methods described in this specification. After the linear equations have been solved, there is a test in block 58 to determine if all non-linear terms have converged in the finite difference equations. If not, the simulator returns to block 54. If the non-linear terms in the finite difference equations have converged, the simulator moves to block 60 and updates values to complete the step in current time. In block 62, the simulator tests to determine if the desired end time (ie stop time) has been reached in the simulation. If not, the simulator advances in time to the next level block 64 and then returns to block 52 to read new variation data over time and start the next step in time. The end point of the simulator has been reached, then the simulator completes the output operations and the operation is determined, block 66. With reference to Figure 10, it must be remembered from our previous discussion that "flow equations" in tank simulators In general, they are discretized using the method of "Finite Volume" (or finite difference) but "Effort Equations" are generally discretized using the "Finite Element" approach. As a result, in Figure 10, the eclipse 46b simulator program includes the following program block: a "Discretization Method for Flow Equations" 68. Furthermore, recalling from our previous discussion that "it would be natural to ask whether finite-volume methods they can be derived for the equations of effort "; and as a result, a "finite volume discretization of a set of stress equations" is presented in this specification. Agree with this, in Figure 10, the eclipse simulator program 46b Figures 5, 6 and 7 includes the following additional program block: a "Finite Volume Method for Linear Elasticity Equations" 70 (ie, "Effort Equations"). In Figure 10, the stage of "Finite Volume Method for Linear Elasticity Equations" 70 receives and responds to the stage of "Method of discretization for Flow Equations" 68. Since the "Flow Equations", in simulators of Commercial deposits are generally discretized using the finite difference or the "Finite Volume" method and the "Effort Equations" are generally discretized using the "Finite Element" approach in Figure 10, the "Finite Volume for Equations" method. of linear elasticity "70 described in this specification represents a possible" Method of discretization for stress equations "using the method of" Finite Volume ". That is, the "Finite Volume Method for Linear Elasticity Equations" 70 of Figure 10 represents a possible implementation of a "finite volume method for linear elasticity" (involving "displacements" and "efforts"), indicating its convenience in meshes of points corners and unstructured in general and demonstrating their order of precision. In Figure 10, the "Finite Volume Method for Linear Elasticity Equations "70 represents a possible" Method of discretization for stress equations "using the method of" Finite Volume "including the following two stages: (1) Interpolation of displacements from cell centers to grid nodes, stage 70a , and (2) Integration over a control volume, to form the discrete system of equations, step 70b, each of the two aforementioned stages 70a and 70b illustrated in Figure 10, associated with the "Finite Volume Method for Equations". Linear Elasticity "70 of Figure 10 described in this specification will be discussed in detail below with reference to Figures 10 to 23. The specification describes a" finite volume approach for discretization of the linear elasticity equations in a general unstructured grid in three dimensions. "This arises from the need to modulate fields of effort in combination with half-pore flow. or multi-phase in a commercial warehouse simulator. In this context, it is naturally convenient to solve the stress equations in the same grid as the one used for the fluid flow and preferably with rock displacement variables located in the same points as the flow variables, such as pressure and saturation. . The grids commonly used in deposit simulators, they are already based on point-corner geometry or are totally unstructured. In any case, they are capable of representing characteristics such as faults, narrowings and refinements of the local grid. In practice, they are often also highly non-orthogonal. It is therefore important to design a discretization of the stress equations that can be implemented in grids of this type. The linear elasticity equations in general domains, are more commonly solved by the finite element method, with displacements located at the vertices of each cell of the grid. However, the flow evolution in a deposit simulator is rarely modeled with finite elements, but in general it uses a finite volume, or finite difference approach. It is natural to ask, therefore, whether the coupled equations can be modeled satisfactorily using a common volume discretization approach for both systems. In Figure 10, this specification provides a partial answer to this question by describing in detail any possible implementation of the "Finite Volume Method for Linear Elasticity Equations" (step 70) of Figure 10, indicating its convenience in corner knit meshes unstructured in general and demonstrating its precision order. Finite volume methods have long been in use in computational fluid dynamics, for discretization of conservation law systems such as the Euler and Navier-Stokes equations (see references 4, 5, and 6 below). The deposit flow equations have also been solved using finite volume methods (see references 7 and 8 below). Applications of finite volume methods to stress equations are rare, but include Demirdzic and Muzaferija (1994) and Demirdzic and Martinovic (1993-see references 9 and 10 below). The last one is aimed at thermo-elasto-plastic stress analysis. Equilibrium stress and reservoir flow equations have been solved within the context of a commercial deposit simulator by Stone et al (2000), Stone et al (2003), and Liu et al (2004-see references 1, 2, and 4 following), which describe a finite element approach to stress equations. Linear Elasticity Equations The equations that describe the linear elasticity of rocks in a hydrocarbon deposit are ? ·? + ?. > = () (i) where b is a vector of body forces and the tensor of effort is given by Here, "u", "v", and "w" are the displacements of unknown rocks, while "p" and "T" are the pressure and temperature resolved during the flow simulation of the deposit. The constants Lamé and, define properties of the rock and can be derived from the Young's modulus and the Poisson's ratio. The Biot constant and the coefficient of thermal expansion T, provide the coupling between the deposit flow and the rock moment equations. If temperature effects are not included in the deposit simulation, then T can be taken as zero. This system of partial differential equations is elliptical, and requires three boundary conditions to be given at all points on the boundary.

Two types of boundary conditions are known: (1) prescribed displacements - the displacements u, v and w are given in a section of the boundary; and (2) prescribed tractions - the traction vector -n is given in a section of the boundary where n is a normal unitary outward direction vector. However, to ensure the existence of a single solution, there are certain restrictions on the type of boundary conditions that may apply. This is due to the fact that the prescribed traction condition defines only displacement derivatives. Clearly, if tractions are prescribed across borders, an arbitrary constant can be added to the offsets for any valid solution and the result will remain as a solution. This can be rectified by insisting that at least one point on the border has a displacement boundary condition. In fact, the selection of border conditions is slightly more restricted than this. At least three different points on the border must have specific displacements and these points must not be collinear. Finite Volume Method In Figure 10, this section of the specification will describe in detail a new "Finite Volume Method for Linear Elasticity Equations" (step 70 of Figure 10) described above. He The method involves two distinct stages: (1) Interpolation of displacements from cell centers to grid nodes, step 70a of Figure 10; Y (2) Integration on a control volume, to form the system of discrete equations, step 70b of Figure 10. These steps are described separately in the following two sub-sections. This is followed by a discussion of the implementation of boundary conditions and coupling of the flow equations. For clarity, most of the description will be for the two-dimensional method, but a section is included that describes some aspects specific to its implementation in three dimensions. Finally, the properties of the method are summarized. Although this request is directed to linear elasticityIt is important to note that the approach described here also applies to many other problems, including the deposit flow equations themselves. In this sense, it can be considered as an alternative to higher order schemes, such as the multi-point flow approach. Interpolation to Grid Nodes, Step 70a of Fig. 10 Referring now to Fig. 11. In Fig. 11, consider a two-dimensional grid section as shown in Fig. 11.

At the center of each cell, such as point 1, an approximation of the two displacements "u" and "v" is stored. As part of the finite volume schema, it is necessary to define a method of interpolating shifts from cell centers to grid nodes. We select to define an interpolation that is accurate for linear functions in any grid geometry. This will ensure that the entire finite volume schema also has this property. In Figure 11, to interpolate grid node A, we will use cell center values in cells 1, 2, 5 and 9 that contain A as a vertex. It will be possible to adjust a polynomial of a higher order than the linear one to these four points, since a linear function in two dimensions only has three free parameters. However, to allow more general unstructured grid geometries, it is convenient to select a linear approximate, and find the best least squares linear fit at the four cell center values. Now with reference to Figures 12 and 13. For example, the grid sections in Figures 12 and 13 contain narrowing cells and local grid refinements show that the number of neighboring cells is variable. Defines the linear function as the form of some function f in the neighborhood of A. Interpolation of all cell center values of neighbors n will require solution of the system of n equations in 3 unknowns.

The least squares solution system is obtained by solving the system equations in 3 unknowns.

ArAe =? T? (5) The value of a obtained from vector c then gives the interpolated value of at point A. If the underlying function was undoubtedly linear, then all n equations are satisfied by the least-squares solution. Therefore, the interpolation is accurate for linear functions in any grid. For grid nodes that are in the border, interpolation may be required to be modified. If displacement boundary conditions are present, they include in the least squares system described above. If the traction boundary conditions occur, there may be less than three neighboring cell center values, which are not sufficient to define a linear interpolant. The solution is direct. Additional adjacent cell center values from the interior are added into the system, until there are at least three non-collinear points. The flexibility of the least squares approach, which can handle an arbitrary value of n > 3 makes this easy to implement. Integration Over a Control Volume, Step 70b of Figure 10 In Figure 11, having obtained the weights required to interpolate from cell centers to grid nodes, the finite volume method advances by integration of the partial differential equations on each grid cell. This process is described for cell O in Figure 11, with vertices ABCD and boundary ^ 'O.

Integrating the equations on O and applying the divergence theorem gives j¡ (v? + b) dV = j? ndS + jfbdV = 0 (6) n 8Ci n which can be expressed as: Approximating the line integrals using the midpoint rule leads to: + FQdyca-GQdxCB (8) * FsdyAD -GsdxAD + Kn, = 0 where Fp denotes F evaluated on the middle side P, dxBA = xB xñ, dyBA = yg and A, etc., and ^ n is the volume of O. It should be noted that the approximation of line integrals by the midpoint rule will be exact if the terms in the stress tensor are linear.

Terms such as Fp involve derivatives of the displacements u and v in the middle-side point P. In order to evaluate these, we configure another linear least-squares interpolant that involves displacements at the two ends of the front and the two adjacent cell-secti-centers. It should be noted that a cell face in any unstructured grid in two dimensions always has two ends and two adjacent cells. The method can therefore accommodate all the required mesh characteristics, such as faults, narrowings and local grid refinements. The least squares interpolant is configured in the same way as for nodal interpolators. Reference is now made to Figure 14. Figure 14 shows the side AB of the control volume ABCD of Figure 11. The interpolant at point P is defined in terms of displacements at points A, B, 1 and 2. Offsets in A and B by themselves are defined as interpolating values from neighboring cell centers as previously described. Again we define a linear function that now focuses on P. The system of interpolation equations is now It should be noted that, unlike the previously described interpolant, there is now a fixed number of four equations in three unknowns. However, this system is resolved in the same way as before. The derivative x of? in P is then obtained in terms of? 1,? 2, ?? and ??, taking the values of q that are obtained from this system. Replacing the expansions of ?? Y ? B in terms of its values of neighboring cell centers, then gives an expression where? x in P in terms of neighboring cell center values. The derivatives and are similarly defined using the value of r. These derived expressions are then used to construct the Fp and Gp flows in P. Boundary Conditions If displacement boundary conditions prescribed in P exist, the boundary condition equation is added to the least squares system to interpolate offsets to grid nodes. If predetermined traction boundary conditions are given, then there is no need to calculate the derivatives of the displacements, as described above. The determined traction vector is simply used directly in the equation for the control volume. An alternate approach to include boundary conditions of displacement, only at the control volume level, and not during interpolation to grid node, was also considered. This makes cell node interpolation more flexible, since it can be used for different amounts of displacements, but it was found to give less precise solutions and produce linear systems that were more difficult to solve. Coupling to the Flow Equations As discussed above, the stress tensor for the voltage-coupled or coupled equations also includes terms that involve pressure and temperature, as well as the Lamé, Biot, and thermal expansion coefficients. These amounts are not considered to be located in cell centers, as is usual in deposit simulation. The treatment of these terms in the current discretization is direct. Each required quantity is simply interpolated to the middle-face point P. The terms can then be included directly in the stress tensor and therefore in the contribution to the equation of the control volume for that face. Implementation in a Regular Grid It is instructive to consider the implementation of this method in a grid regulating the step length h. For simplicity, we continue to consider only the two-dimensional case and restrict attention to pure Linear Elasticity without coupled flow terms. The conclusions lead to the most general case. In a regular grid, the moment equation x has the template or stencil: If we consider this stencil applied to an infinite grid in the absence of boundary conditions, it is clear that the constant displacements u and v satisfy the homogeneous equation with zero force of the body. However, there are no spurious trivial solutions. This property, known as discrete ellipticity, has important consequences. It means that the solution of the discrete system in a genuine mesh will not be prone to non-physical oscillations. Besides, the Discrete ellipticity also implies that the solution of the linear system by an iterative method such as multiple grids, in general will be efficient. There are alternate finite volume discretizations that do not have this important property. One obvious is to use the trapezoidal rule instead of the midpoint rule to calculate the flow on the AB side of the ABCD control volume: *. { + FA ½AD ~ i. { D + GA) dxAD 4- Va This produces the stencil This stencil still admits the trivial constant solution, but also admits an oscillating solution of spurious chessboard. For this reason, it is not discretely elliptical and will produce non-physical oscillatory solutions and it will be difficult to solve. Three-dimensional implementation The implementation of the finite volume method in three dimensions, follow the same lines described above for two dimensions. Interpolation to grid nodes naturally uses a three-dimensional linear function, with four undetermined coefficients instead of three. This requires at least four neighboring cell-center values and results in a 4-by-4 linear-squared minimum-square system solution. Near borders, additional interior points are included as for the two-dimensional case. The integration of control volume is rather complicated, since the faces between neighboring cells in three dimensions are polygons in a three-dimensional space instead of simple lines. These polygons in general are not in a single plane. The integration of the flow on these polygons is handled by defining a polygon center, and dividing the polygon into a set of triangles defined by the center and each pair of connected polygon vertices. The flow over each triangle is then added to define the flow over the entire face. Properties of Discretization The finite volume discretization of the voltage-coupled or coupled equations described above has the following properties. (1) When applied to a regular grid, the discretization has a compact 9-dot stencil in two dimensions and a stencil of 27 points in three dimensions. Nearby borders with traction border conditions, this stencil can be locally extended, but easily accommodated in a sparse matrix storage format. (2) Unknown displacements are stored in cell centers; as well as the flow variables in the simulation of the coupled deposit. (3) The interpolation to grid nodes is accurate for linear functions, while the integration of control volume is accurate for a linear stress tensor. The total discretization is therefore accurate for linear displacements and a linear stress tensor. (4) The method is conservative by construction, as is generally the case for finite volume methods. That means that the flow on one side of the grid is the same when it is calculated on both cells that contain it. Consequently, if the equations are integrated numerically over the entire domain, the Internal flow contributions are canceled and the resulting integral depends only on border flows. This is the discrete equivalent of the divergence theorem used to construct finite volume discretization. (5) The method can handle general unstructured grids, including mesh characteristics such as faults, narrowings and local grid refinements, (6) The method is discretely elliptical. Numerical Results In order to validate that the current method is accurate for linear solutions in any grid and to establish its precision order, it was implemented in general unstructured meshes and tested for a number of known real solutions in a sequence of refined meshes and meshes with particular characterisitics. Validation of Linear Solutions Accuracy. In this case, the real solution was chosen to be the set of linear functions: u = 3x ~ y + 2z + l v = -2x-4y + z + 6 (13) w = x + y + z -3 The pressure was also chosen to be linear, with a constant temperature: /? = 1000 * + 2000.y - »- 3000z + 4000, T = 4000 (14) The various coefficients were: Á = x + y- z + 3, μ = 3? + 4? + 10, a = 1, at = 1 (15) Since the actual stress tensor is linear for this solution, the finite volume method must be accurate in any mesh. This was tested for the following cases: (1) A sequence of regular three-dimensional meshes with increasing refinement. (2) The same sequence with each mesh node randomly disturbed. (3) The same sequence mapped by coordinate transformation in a ring. (4) Meshes with faults, narrowings and local grid refinements. (5) Meshes transformed with faults, narrowings and local grid refinements. Reference is now made to Figures 15, 16, 17 and 18. In Figures 15-18, in all cases, the Error rule between the numerical solution and the actual solution was calculated and found to be effectively zero. Figures 15, 16, 17 and 18 illustrate the linear nature of the displacements in sliced planes, through the various geometries tested. Convergence order for known general solution In this case, a known set of displacements is chosen arbitrarily as: U -senp x) cos (r y / 2) sin (Vr z) v = 5x2 + 2y2 + 3z2 + \ (16) w = elz cos (+ y) The pressure, temperature and various coefficients are all chosen as for the linear test case. For this general solution, the stress tensor is not linear and the body forces, therefore, are not zero. The actual values of these terms were calculated and included on the right-hand side of the discretization. Reference is now made to Figure 19. The problem was solved in the same mesh sequence described above. The accuracy of the numerical solution was measured when calculating its error standard. Figure 19 shows the logarithm of the error standard in a mesh sequence, for various transformations of mesh. The slope of the graph gives the order of convergence of the finite volume method, since the graph is the logarithm of the error standard against the logarithm of the step length. More negative values of log (h) correspond to finer meshes. In this region, Figure 19 clearly indicates a slope of approximately 2 for all three geometries, i.e. a cube, a cutting region and a three-dimensional ring. This indicates that the infinite volume method is of second order precision in quite general meshes. The dispersed system of linear equations is solved using the SAMG algebraic multiple grid method (see reference 11). Since the system is discretely elliptical, we expect this to be an efficient solution technique and to demonstrate independent mesh convergence speeds. Reference is now made to Figure 20 illustrating "Table 1". In Figure 20, "Table 1" shows the number of multi-grid interactions (with upe times in parentheses) required to solve the system in a sequence of meshes with N nodes in each of three dimensions. The number of equations in the system linear is shown as neqn. As can be seen in Table 1 in Figure 20, the number of interactions of multiple grids required is quite independent of both the number of equations and the mesh transformation. This is possible optimal convergence behavior and is characteristic of multiple grids applied to a discretely elliptical problem. Similar experiments in the previously mentioned alternating non-discrete elliptical method showed a significant increase in the number of iterations required as mesh was refined. Integration within a commercial simulator The finite volume method discussed in this specification has been implemented within a commercial warehouse simulator, ie the Eclipse simulator 46b shown in Figures 5, 6 and 7. A partially coupled mesh-in-motion approach it is used where the flow equations are advanced to a reporting stage and the resulting pore temperature and pressure are used as source terms in the stress tensor for calculating the subsequent stress field. The resulting stress is used to derive an updated porosity using the porosity-spherical relationship and the flow simulation advances to the next stage of report. The same coupling approach has been used for the geomechanical option of existing finite elements (see Reference 3 below). Reference is now made to Figures 21, 22 and 23. In Figures 21 and 22, results for finite volume and finite element discretizations for a simple test example running within the simulator are illustrated in Figure 21 and the Figure 22. These show good qualitative and quantitative agreement. One difference is in the formation of stress layers near the wells (which do not complete through the depth) in the case of finite volume. This is under investigation. Figure 23 shows a case of infinite volume that includes local grid refinements in the vicinity of the wells. Comparing the two approaches, the main differences and similarities are as follows: (1) The finite volume method has somewhat smaller independent variables, which are located with the flow variables in grid block centers, instead of grid nodes. This does not necessarily imply faster solution times, although it can simplify the process of coupling the flow and effort equations more implicitly. (2) The finite volume method aims to manage all the mesh characteristics generically through its unstructured mesh formulation. The finite element approach currently implements restrictions and refinements of the local grid that are in the implementation process. Faults have not yet been included. Nevertheless, there is a general aspect in the generic nature of the finite volume approach. The cost of configuring the grid connectivities at the beginning of the simulation generally makes the initialization of the finite volume method slower than the finite element counterpart. (3) The finite volume method is conservative. This means that the flows in any given inner face of the grid are calculated exactly in the same way for each of the two cells that contain it. If the equations are integrated numerically over the whole domain, therefore the internal flows are canceled, and the result depends only on the border flows. This is a discrete analog of the divergence theorem. (4) Inactive cells are not included in the infinite volume model but are included as fictitious or test equations for finite elements. (5) Both methods produce large systems scarce linear ones, which are solved by the solver or mathematical mechanism of solution of multiple SAMG grids. The resolution times required are similar in both cases. (6) Both methods retain good accuracy in general non-orthogonal meshes. The previous finite difference approach1 does not have this important advantage. A "new Finite Volume Method for Tank Flow Equations and effort coupled in three dimensions" has been presented. This method exhibits many desirable characteristics such as second order accuracy in general unstructured meshes including faulted meshes, restrictions and refinements of local grid. The method has been shown to be discretely elliptical and therefore to produce solutions without non-physical oscillations and be suitable for efficient solution using multiple grid methods. The convergence properties have been demonstrated by application to a sequence of model test problems with known solutions. Nomenclature a = Biot constant, dimensionless. aT = coefficient of linear expansion, 1 / degrees C. , 1 / degrees F? = Lamé constant, bar, psi, Pa ss = normal voltage, bar, psi, Pa μ - Lamé constant (modulus of stiffness), bar, psi, Pa REFERENCES The following References include References 1 to 11 set forth below, each is incorporated by reference into this specification: 1. Stone, T., Bowen, G., Papanastasiou, P. and Cook, J., "Coupled geomechanics in a commercial reservoir simulator", SPE 65107 presented at SPE EUROPEC 2000 Petroleum Conference, October 2000. 2. Stone, T., Xian, C, Fang, Z., Manalac, E., Marsden, R., Fuller, J., "Coupled geomechanical simulation of stress dependent reservoirs ", SPE 79697, filed at SPE Reservoir Simulation Symposium, February 2003. 3. Liu, Q., Stone, T., Han, G., Marsden, R. , Shaw, G., "Coupled stress and fluid flow using a finite element method in a commercial reservoir simulator", SPE 88616, presented at SPE Asia Pacific Oil and Gas Conference in Perth, Australia, October 2004. 4. Peyret, R, and Taylor, TD, Computational methods for fluid flow, Springer Verlag, Berlin, 1983.

. Jameson, A., Schmidt, W. and Turkel, E., "Numerical solutions of the Euler equations by finite volume methods using unge-Kutta time stepping", AIAA Paper 81-1259, 1981. 6. Crumpton, PI and Shaw, GJ, "Cell vertex finite volume discretizations in three dimensions", International. Journal for Numerical Methods in Fluids, vol 14, 505-527, 1992. 7. Rozon, B. J., "A generalized finite volume discretization method for reservoir simulation", document SPE 18414 presented at the Reservoir Simulation Symposium in Houston, Tex. , feb. 6-8, 1989. 8. Crumpton, P. I., Shaw, G. J. and Ware, A. F., "Discretization and multigrid solution of elliptic equations with mixed derivative terms and strongly discontinuous coefficients", Journal of Computational. Physics, vol 116, 343-358, 1995. 9. Demirdzic, I. and Muzaferija, S., "Finite volume method for stress analysis in complex domains", International Journal for Numerical Methods in Engineering, vol 37, 3751-3766, 1994. 10. Demirdzic, I. and Martinovic, D., "Finite volume method for thermo-elasto-plastic stress analysis", Computer Methods in Applied Mechanics and Engineering, vol 109, 331-349, 1993. 11. Stuben, K., SAMG User 1 s Manual Reléase 21bl, Fraunhofer Institute SCAI, Schloss Birlinghoven, D-53754 St Augustin, Germany, July 2002. The above description of the "Finite Volume Method for Linear Elasticity" being thus described, will be evident that it can be varied in many ways. Said variations shall not be construed as a separation of the spirit and scope of the method or apparatus or storage device of claimed programs, and all such modifications as will be apparent to a person skilled in the art, are intended to be included within the scope of the following claims.

Claims (3)

1. CLAIMS 1. A method implemented by computer to perform an effort calculation adapted to model a set of efforts and displacements in a deposit, characterized in that it comprises: (a) building a model of a deposit on a region of interest based on a plurality of the records of exit of well logs and seismic data by applying a workstation processor to form grid of the region of interest forming a grid, the grid is constituted by one, cells and has nodes, each cell has a center of cells where the grid stores in a storage medium; (b) interpolate unknown rock displacements in the region of interest from cell centers to grid nodes; (c) integrate over each cell to form a discrete system of equations relating to the deposit, where the described system of equations is stored in the storage medium; and (d) use the system of discrete stored equations to model efforts and displacements in the deposit.
2. 2. The computer implemented method according to claim 1, characterized in that the step of interpolation (b) to interpolate unknown rock displacements in the region of interest of cell centers for grid nodes comprises: (bl) defining a linear function ^ = a + b (x xA) + c (y yA) as a form of a function f in a neighborhood of A, where A is a grid node; (b2) interpolate all neighboring cell center values. "3. The computer implemented method according to claim 2, the step of interpolation (b2) comprises: solving a system" with three unknowns as follows: 4. The method implemented by computer according to claim 3, characterized in that a solution of least squares of the system of n equations in three unknowns is obtained by solving a system of three equations with three unknowns as follows, ATAc = AT, where c is a vector, and where a value of a obtained from the vector c produces an interpolated value of f in grid nodes A. 5. The method implemented by computer in according to claim 1, characterized in that the integration stage (c) to integrate on each cell to form a discrete system of equations, comprises: (the) integrate the system of "equations with three unknowns on a cell O and apply a theorem of divergence producing in this way, JJ (V -? + B) F = j? - a dS + jb dV = 0 O. O that can be expressed as where 6. The computer implemented method according to claim 5, the integration step (c) further comprises: (c2) approximating a set of line integrals using a midpoint rule producing this way, * FQdyCB ~ GQdxCR ± F * dyDC-GRdxDC where Fp denotes F evaluated in point of the middle side P, is a volume of cell O. 7. A machine-readable program storage device that tangibly incorporates a set of instructions executable by the machine to perform steps of the method to conduct an effort calculation adapted to model a set of efforts and displacements in a deposit; method steps are characterized in that they comprise: (a) instructions for constructing a deposit model of a deposit on a region of interest based on a plurality of logs and seismic data output records upon application of a workstation processor to form grid of the region of interest, the grid comprises one or more cells that have nodes, each cell has a cell center where the grid is stored in a storage medium; (b) instructions for interpolating unknown rock displacements in the region of interest from cell centers to grid nodes; (c) instructions to integrate on each cell to form a system of discrete equations concerning the deposit, where the system of discrete equations is stored in the storage medium; and (d) instructions for using the discrete stored equation system to model stresses and displacement in the reservoir. The program storage device according to claim 7, characterized in that the step of interpolating (b) to interpolate unknown rock offsets in the region of interest of cell centers for grid nodes, comprises: (bl) instructions for define a linear function ^ = a + b (x xA) + c (y yA) as a form of a function f in a neighborhood of A, where A is a grid node; and (b2) instructions to interpolate all n neighboring trap center values. 9. The program storage device according to claim 8, characterized in that the step of interpolation (b2) comprises: solving a system of "equations with three unknowns as follows: 10. The method implemented by computer according to claim 9, characterized in that a least squares solution of the system of n equations in three unknowns is obtained by solving a system of three equations with three unknowns as follows, ATAc = AT #? , where c is a vector, and where a value of a obtained from vector c produces an interpolated value of f in a grid node A. 11. The method implemented by computer in accordance with claim 7, characterized in that the step of integration (c) to integrate on each cell to form a system of discrete equations, comprises: (the) instructions to integrate the system of "equations with three unknowns on a cell O and apply a divergence theorem producing in this way, j] (V -? + b) 7 - c [? . n < S + j¡bdV = 0 that can be expressed as where 12. The storage device of the program according to claim 11, characterized in that the integration step (c) further comprises: (c2) instructions for approximating a set of line integrals using a mid-point rule producing this way, -GQdxcg where Fp denotes F evaluated in P medium side, dxBñ = xB xA, dyBA = yB yA and? O is a volume of cell O. 13. A system adopted for conducting an effort calculation for modeling a set of efforts and displacements in a tank, characterized in that it comprises: a first apparatus adapted to construct a deposit model of a deposit on a region of interest, based on a plurality of well log output record and reduced seismic data by applying a workstation processor to form grid of the region of interest, the grid comprises one or more cells and has nodes, each cell has a cell center where the grid is stored in a storage medium; a second apparatus adapted to interpolate unknown rock shifts in the region of interest from cell centers to grid nodes; a third apparatus adapted to integrate on each cell to form a system of discrete equations relating to the deposit, wherein the discrete system of equations is stored in the storage medium; and a fourth apparatus adapted to model the reservoir while using the system of discrete stored equations to model stresses and displacement in the reservoir. The system according to claim 13, characterized in that the second apparatus is adapted to interpolate unknown rock displacements in the region of interest from cell centers to grid nodes, comprising: an apparatus adapted to define a linear function ^ = a + b (x xA) + c (y yA) as a form of the function f in a neighborhood of A, where A is a grid node; and the apparatus is adapted to interpolate all neighboring cell center values 15. The system according to claim 14, characterized in that the second for problems comprises: an apparatus equipped to solve a system of n with three unknowns as follows: 16. The system according to claim 15, characterized in that a least squares solution of the system of n equations with three unknowns is obtained by solving a system of three equations with three unknowns, as follows, ATAc = AT < p, where c is a vector, and where a value of a obtained from the vector c produces an interpolated value of f in the grid node A. 17. The system according to claim 13, characterized in that the third apparatus is adapted to integrate over each cell to form a system of discrete equations, comprising: an apparatus adopted to integrate the system of "equations with three unknowns on a cell O and apply a divergence theorem producing in this way, which can be expressed as where F - 18. The system according to claim 17, characterized in that the third apparatus further comprises: an apparatus adapted to approximate a set of line integrals using a mid-point rule producing this way, FptyaA - G? DxM + FsdyAO -? ?? ^ + Vab t = 0 where Fp denotes F evaluated in the P of the middle side, dxM - xB ~ xA ^ dyQÁ ^ yB -yAj and V is a volume of cell O. 19. A computer program adapted to be executed by a processor, the computer program when executed in the processor practices a process adopted to conduct a calculation of effort to modulate a set of efforts and displacements in a deposit, the process is characterized because it comprises: (a) constructing a deposit model of a deposit on a region of interest based on a plurality of records of output from well logs and seismic data by applying a workstation processor to grid the region of interest, the grid comprises one or more cells and has nodes, each cell has a cell center where the grid is stored in a storage medium; (b) interpolate unknown rock displacements in the region of interest from cell centers to grid nodes; (c) integrate over each cell to form a system of discrete equations concerning the deposit, where the discrete system of equations is stored in the storage medium; and (d) use the system of discrete stored equations to model efforts and displacements in the deposit. The computer program according to claim 19, characterized in that the step of interpolating (b) to interpolate unknown rock displacements in the region of interest of cell centers for grid nodes, comprises: (bl) defining a linear function ^ = a + b (x xA) + c (y yA) as a form of a function f in a neighborhood of A, where A is a grid node; and (b2) interpolating all neighboring trap center values 21. The computer program according to claim 20, characterized in that the step of interpolation (b2) comprises: solving a system of n equations with three unknowns as follows: 22. The computer program according to claim 21, characterized in that wherein A solution of minimums One square of the system of "equations in three unknowns is obtained by solving a system of three equations with three unknowns, as follows, ATAc = AT <p, where c is a vector, and where a value of a obtained from the vector c produces an interpolated value of f in the grid node A. computer program according to claim 19, characterized in that the integration step (c) to integrate on each cell to form a system of discrete equations, comprises: (the) integrate the system of n equations with three unknowns on a cell O and apply a divergence theorem producing in this way which can be expressed as where 24. The computer program according to claim 23, characterized in that the integration step (c) further comprises: (c2) approximating a set of line integrals using a mid-point rule thus producing ^ Qdyca-GQ xCB + FsdyAD - + Vb = 0 where Fp denotes F evaluated in P of the middle side, is a volume of cell O. 25. The computer program according to claim 1, characterized in that the step of using the discrete system of equations to model the forces and displacements in the deposit comprises displaying simulation results to a viewer, based on a system of equations discreetly stored. 26. The computer program according to claim 7, characterized in that the step of using the discrete system of equations for modeling the efforts and displacements in the deposit, comprises displaying simulation results to a viewer based on the system of discrete stored equations. 27. The system according to claim 13, characterized in that the fourth apparatus is adopted to use the discrete system of equations to model the forces and displacements in the tank., comprises an apparatus adapted to display simulation results to a viewer, based on the system of discrete stored equations. 28. The computer program according to claim 19, characterized in that the step of using the discrete equation system to model the forces and displacements in the deposit comprises displaying the simulation results to a viewer, based on the system of discrete stored equations.