NO344113B1 - Efficient application of reduced variable transformation and conditional stability testing - Google Patents

Description

Description

The subject area of the invention

The invention generally relates to computer-activated reservoir simulation of fluid flow in the underground reservoir and in particular to composite reservoir simulation.

Background

Conceptually, complex reservoir simulation to simulate flow in an underground, hydrocarbon-bearing reservoir can be considered as a modeling of a number of interconnected mixed tanks of fluids (cells) at a given pressure, a certain temperature and compositions. As time passes (the simulator can time step towards a final point in time and the results are searched), the conditions in the tanks change as a result of the fluid movement, wells and other external factors. Flash calculations are necessary to determine for each new set of pressure, temperature and overall fluid composition, the number of fluid phases and their sizes and compositions. This calculation fundamentally involves finding the smallest thermodynamic state function (Gibbs Free Energy (GFE)) which is periodic in nature and consequently difficult to collect and computationally expensive, especially when detailed fluid models are used, i.e. when there are many hydrocarbon components. It is therefore of great interest to be able to develop algorithms that are computationally efficient and yet reliable and accurate.

For the purpose of the description, the activity carried out in "flash" calculation shall be further divided into stability testing which can reveal instability in a given phase of the current conditions and "split" calculations which aim to determine balanced phases and compositions for an assumed phase configuration.

'describes a method for increasing the efficiency of flash calculations on, in a computational reservoir simulator, is described by Claus P. Rasmussen, Kristian Krejbjerg, Michael L. Michelsen and Kersti E. Bjurstrom, Increasing the Computational Speed of Flash Calculations with Applications for Compositional , Transient Simulations, Society of Petroleum Engineers, SPE 84181, February 2006 SPE Reservoir Evaluation & Engineering. While performing flash calculations, most of the time is spent performing stability analysis. Rasmussen et al. proposed a criterion to skip many of the stability analysis checks.

'Three-Phase Flash Calculation for Hydrocarbon Systems Containing Water' by S. Mokhatab (2002) describes a flash algorithm which is based on the thermodynamic condition of equal fugacies for each component in each phase. The resulting phase compositions then provide better values for updating the distribution coefficients, using an equation of state. Equations of state (EOS) have been successfully used to describe the phase behavior of reservoir crude and gas condensates, but the phase behavior of water/reservoir crude oil has not yet been predicted by an EOS. In this study, it is shown that using the current algorithm with the Shinta and Firoozabadi association model gives reliable and better results in comparison with experimental data for three-phase calculation in the VLLE region In fig. 1 shows a pressure/temperature map for a fluid in a cell. Point A is shown in the two-phase area where both a gas phase and a liquid phase exist. Point B is located on a transition line between a two-phase area and a one-phase area (the phase boundary). Point C is located in a "shadow zone" of single-phase areas close to the two-phase area. Finally, point D lies in a remote area of a phase domain. A vertical line is also shown separating the single-phase liquid on the left and the single-phase gas on the right. Depending on where the estimated phase state of the fluid in a cell is located, some stability calculations can be omitted instead of performing a stability analysis for all cells during the iterations in a time period. This general criterion for converting calculations shall be described in detail below in connection with stability testing in Part 5.

The sub-steps corresponding to the "split" and "stability" calculation described in Rasmussen et al., use a traditional method, which used the solution of non-linear problems with a size corresponding to the number of hydrocarbon components. There is scope to improve the efficiency of these sub-steps, especially for simulation models that involve a large number of components.

Firoozabadi, A. and Pan, H., Fast and Robust Algorithm for Compositional Modeling: Part I – Stability Analysis, SPE 63083 and Firoozabadi, A. and Pan, H., Fast and Robust Algorithm for Compositional Modeling: Part II – Two- Phase Flash, SPE 71603, describes the application of reduced variable strategies for stability and split calculations in a complex reservoir simulation, but the authors do not describe how the stability tests can be avoided. Also, certain stability algorithms that are formulated may experience convergence difficulties, especially under conditions that extend into the undersaturated zone. In addition, the splitting algorithm is formulated for the vapor phase and will experience numerous and/or convergence difficulties near the dew points due to the virtually non-existent liquid phase.

Newton's method is usually used to solve non-linear systems of equations. However, it must be ensured that the repetitions do not exceed the physical limits of unknowns. By using Newton's method for the problems of phase behavior formulated in the form of reduced variables, there is a need to ensure that physical limits on the reduced variables are not violated.

The disadvantages of previous methods for composite reservoir simulations mentioned above will be discussed below.

A method, system, and computer-readable media that executes instructions for performing a complex reservoir simulation of a subsurface hydrocarbon-bearing reservoir is provided. Reduced variables for flash calculations are utilized together with a method of state stability testing to be able to achieve an optimal efficiency of phase behavior calculations in a composite reservoir simulator. Preferably, a least abundant phase is also selected as primary variables associated with a primary phase and a secondary phase are selected for a more abundant phase, thus ensuring stability by not dividing by a value close to zero due to the selection of the primary phase becomes associated with a phase that is not the least abundant. A limit interval can also be used to limit solution changes in reduced, variable algorithms (phase split and stability) in order to achieve a greater stability of the algorithms. Furthermore, stability tests can be performed during flash calculations using reduced variables and indirect existing form based on the definition of reduced variables and the tangent plane distance condition.

It is an object of the invention to combine the concept of reduced variables for flash calculations with a method of condition stability testing in order to achieve optimal efficiency of phase behavior calculations in a composite reservoir simulator.

It is the purpose of the invention to provide a more reliable, reduced variable phase splitting algorithm by selecting main variables corresponding to the least abundant phase that exists.

It is another object of the invention to use a limit interval to limit solution changes in the reduced, variable algorithms (phase split and stability) in order to achieve greater stability of the algorithms and/or avoid excessive repetitions.

It is further an object to provide an improved method for performing stability tests and using reduced variables and a direct existing form based on the definition with reduced variables and the tangent plane distance condition.

Summary of the invention

In a first aspect, the invention describes a method for composite reservoir simulation of a fluid flow, comprising providing a reservoir model for an underground reservoir where the model has a grid that defines a plurality of cells in such a way that each cell is associated with the fluid therein, where the fluid is in a vapor phase and a liquid phase. The method further comprises selecting at least one cell from the plurality of cells associated with fluid that is in the vapor phase and fluid that is in the liquid phase, estimating which of the vapor phase or the liquid phase is present in the least abundant amount and assigning the phase with the least abundant amount as primary phase and assign the second phase as secondary phase. Furthermore, the method includes calculating phase characteristics of the primary phase which is responsive to primary variables associated with the primary phase and calculating the phase characteristics of the secondary phase using the mass balance and which is responsive to secondary variables associated with the secondary phase, so that stability is ensured during calculation of the phase characteristics of the secondary phase due to allocation of the phase with the least abundant amount as the primary phase. Further, the method comprises emitting a visual display that is responsive to the calculated phase characteristics of the primary phase and to the calculated phase characteristics of the secondary phase.

In a second aspect, the invention describes a method for determining the composition of fluid in a cell during computer-activated reservoir simulation. The method comprises the steps: a) Providing a reservoir model for an underground reservoir where the model has a grid defining a plurality of cells in such a way that each cell is associated with the fluid therein. (b) To perform one of the following calculations for each cell based on the composition of the fluid in the cell at the previous time step: (i) To perform direct reduced variable phase split calculations using K values when the composition of the cell has a fluid of several phases in a previous time step, where the directly reduced, variable phase split calculations are performed so that the phase composition of the fluid that is present in the least abundant amount is the primary phase. (ii) To perform a single-sided, reduced, variable stability test using vapor-initiating modules when the composition of the fluid in the cell comprises a single phase placed in the shadow region, i.e. the liquid side in the previous time step. (iii) To perform a single-sided, reduced, variable stability test using liquid-initiating modules when the composition of the fluid in the cell comprises a single phase located in the shadow region, i.e. the vapor side in the previous time step. (iv) Not performing a stability test when the composition of the fluid in the cell comprises a single phase in the remote region at the previous time step. (c) Determining whether there is an error in step (b) and performing one of the following steps: (i) Performing a first flash calculation to determine fluid composition in the cell when it is determined that there is an error in step (b). (ii) Using the calculated result from step (b) when determined not to be in error in step (b). The method further comprises the step of: (d) Emitting a display representing the composition of the fluid in response to step (c).

In a third aspect, the invention describes a system used during reservoir simulation to determine the composition of fluid in a cell. The system comprises a data receiver which receives input data from a source, where the input data comprises a reservoir model and data. Furthermore, the system comprises a processor and a computer program executable on the processor, the program comprising: (a) A least abundant amount allocation module which (i) calculates whether a vapor phase or a liquid phase of fluid to the cell is present in least abundant amounts in response to input data received by the data receiver and (ii) assigns the phase having the calculated least abundant amount as the primary phase and assigns the other phase as the secondary phase. The program further includes: (b) A fluid phase properties calculator which (i) calculates the phase properties of the primary phase responsive to the primary variables associated with the primary phase, and (ii) calculates the phase properties of the secondary phase using mass balance and secondary variables associated with the secondary phase, the fluid phase properties calculator ensuring stability while calculating the phase properties of the secondary phase due to the least abundant amount, the assignment module assigns the phase with the least abundant amount as the primary phase. Furthermore, the program includes: (c) An output producer adapted to produce and communicate the calculated phase properties of the fluid in a readable format.

In a fourth aspect, the invention comprises a method for determining composition of fluid in a cell during a computer-enabled reservoir simulation, the method comprising the steps of: (a) Providing a reservoir model for an underground reservoir where the model has a grid defining a plurality of cells on such a way that each cell is connected to the fluid therein. (b) Determining whether each cell has a single phase or multiple phases in a previous time step, (i) whether the cell had multiple phases in the previous time step, perform directly reduced variable phase split calculations by choosing K values, (ii) whether the cell had a single phase in the previous step, perform the following steps depending on whether the cell's previous location in a phase plane: (1) if the fluid is in the shadow region, i.e., the liquid side, perform a single-sided reduced variable stability test using vapor-initiating moles, (2) if the fluid is in the shadow region, i.e. the vapor side, perform a single-sided, reduced variable stability test using liquid-initiating moles, (3) if the fluid is in the remote region, do not perform a stability test cell to determine fluid composition and then skip to step (c). (c) To determine whether there is an error in expedited processing at (i), (ii) or (iii). (i) If there is an error, perform an initial flash calculation on the cell to determine fluid composition, (1) if the fluid is found to be a single phase, perform additional calculations to determine location in the phase plane for use in subsequent iterations, (2) if the fluid is found to be in a plurality of phases proceed to step (d). (ii) If there are no errors, proceed to step (d). (d) To emit a display representing the composition of the fluid.

In a fifth aspect, the invention comprises a computer program stored on a computer-readable medium used for composite reservoir simulation of fluid flow in an underground reservoir. The program comprises a least abundant allocation module which estimates which of the fluid in the vapor phase and the fluid in the liquid phase is present in the least abundant amount in a grid cell, and assigns the phase with the least abundant amount to the primary phase and assigns the other phase to the secondary phase . Furthermore, the program includes a calculator for fluid phase properties that calculates phase properties for the primary phase that is responsive to the primary variables associated with the primary phase and calculates phase properties for the secondary phase using mass balance and that is responsive to secondary variables associated with the secondary phase. Furthermore, the program includes an output producer module that communicates the calculated phase properties in a readable format to a visual display.

Brief description of the figures

The invention shall be described in more detail below with reference to the drawings, where

fig. 1 is a pressure-temperature diagram outlining several areas in the phase plane to show the concepts central to condition stability testing,

fig. 2 is a flowchart showing the combined use of conditional stability test logic for the total flash update, and reduced variable algorithms for iteratively solving stability and phase split problems at a particular time step in a composite reservoir simulator,

fig. 3 shows physical limits used for the reduced variables to ensure the Newton repetitions,

fig. 4 is a functional diagram of a non-linear iteration loop with flash calculations made during calculated simulation of fluid flow according to the invention in the context of an underwater hydrocarbon-bearing reservoir model,

fig. 5 is a functional diagram of an embodiment of a method according to the invention,

fig. 6 is a functional block diagram of another embodiment of a method according to the invention, and

fig. 7 is a schematic view of an embodiment of a system and computer readable media according to the invention.

Detailed description

The following nomenclature is used in the following equations:

Symbols

D = Tangent plane distance (TPD)

G = Gibbs Free Energy function (GFE)

c = Number of hydrocarbon components

m = Number of non-zero Eigenvalues

M = Number of reduced parameters. Equal to 1

P = Pressure

Q = Reduced variables (M-vector) of records Qa

Q = Reduction coefficient matrix of size Mc of entries qij

R = Universal gas constant

T = Temperature

x = Liquid phase composition (mole fractions)

y = Vapor phase composition (counter fractions)

Y = Non-normalized moles of sample phase

z = Feed composition (total) (mole fractions)

Subscribe

i,j = Component index

α = Reduced variable index

L = Liquid

V = Steam

Superscripts

F = Feeding phase

T = Trial phase

Greek symbols

Φ<1 = Fugacity coefficient>

δij= Binary interaction coefficients (BIC)

The descriptions from Firoozabadi, A. and Pan, H., Fast and Robust Algorithm for Compositional Modeling: Part I – Stability Analysis, SPE 63083 and Firoozabadi, A. and Pan, H., Fast and Robust Algorithm for Compositional Modeling: Part II – Two-Phase Flash, SPE 7160 is incorporated here in its entirety. Likewise, reference is made to the contents of US patent 2006/0036418 to Pita et al., Highly-Parallel, Implicit Compositional Reservoir Simulator for Multi-Million-Cell Models, in its entirety. Finally, reference is made to Michael Michelsen and Kersti E. Bjurstrom, Increasing the Computational Speed of Flash Calculations with Applications for Compositional, Transient Simulations, Society of Petroleum Engineers, SPE 84181, February 2006 SPE Reservoir Evaluation & Engineering. Finally, the descriptions contained in Michael L. Michelsen, The Isothermal Flash Problem, Part II, Phase-split Calculation, Fluid Phase Equilibria, 9 (1982) 21-40 are also incorporated by reference in their entirety.

Fig. 4 shows the general steps taken during a non-linear iteration loop. Calculations are made about properties and EOS. A Jacobin matrix is then generated. A linear solution method is used to solve a linear set of equations for a solution. The solution is then tested for sufficient convergence. If not sufficiently converged, new calculations are made for EOS and property. Or are the converted results used.

1 Cubic equation about the condition

For a pure substance, an equation of state (EOS) is a mathematical relationship between pressure, temperature and volume, for a mixture a composition is added to this relationship. The third degree form of EOS is the most popular and in particular the Redlich-Kwong-Soave-Peng-Robinson family of EOS has long been the industry standard for complex reservoir simulation.

The preferred EOS is written generically in this print explicit form

(1)

temperature-dependent active expression a = a(T) and where b is the repulsive expression.

The specific EOS corresponding to the parameter choices m1 and m2 is given below:

Defining the compressibility factor, or Z-factor,

<(2)>

GRT

a rearrangement of (1) by using definitions (2)-(4) gives the cubic forms EOS,

which is solved depending on the phase and other conditions for the right root Z =Z(A,B).

By applying EOS to a mixture, as opposed to a pure substance, mixing rules are used to calculate the parameters a and b (or equivalently A and B).

The most commonly used mixing rule for the attractive parameter is the symmetric double sum

Here are the binary interaction coefficients (BIC) which provide for chemical interactions between components of different types. It is a symmetric matrix type with zero diagonal entries.

The supported parameter is usually calculated as a mole average,

In the expressions (6)-(7), the pure component expressions are given by:

which introduces reduced pressure temperatures

designates component critical properties. The functional form of the temperature dependence function depends on the EOS chosen. By using Wifor to name the component of rising factor, the definition becomes the conditions as follows:

The "corrected" form of Peng-Robinson was also supported as an alternative to (12) of the form

The expression Ωbj= Ωboer a constant depending on the temperature for the Redlich-Kwong-Soave-Peng-Robinson family of EOS. The default values of the EOS constants Ωa0 and Ωboer are given in the table below and can be overridden by the user:

A traditional weakness with two-parameter EOS, for example (1) is often poor prediction of the liquid density. To solve this drawback, a three-parameter expansion via the standard Peneloux et al. volume shift can be used. In this connection, the molar volume is calculated in accordance with

Here v<eos>molar volume is predicted by EOS, equations (2) and (5) and the corrected expression is calculated from

where Xj denotes the phase composition and Cier a set of volume changes associated with

user-supplied non-dimensional volume shifts Si according to

Finally, fugacity coefficients and their derivatives are fundamental building blocks in the construction of EOS algorithms. They can be calculated directly from EOS using first principles. For an EOS of the generalized type (1), they can appear in the form

The fucacity of a component is a measure of its "tendency to escape" from a phrase and is therefore directly applicable to equilibrium calculations

The condition for thermodynamic equilibrium can be expressed as

equivalent to φ' from which the definition of K values can be introduced as

ϊ ΐ

2 Reduced variable approximation to EOS

As mentioned above, the matrix is l-dy which is symmetric and consequently provided with a complete, orthogonal set of eigenvectors and corresponding real Eigenvalues. By the spectral expansion theorem, this matrix can therefore be shown as a series

This actual ranking of the matrix 1- δij will depend on the number of non-hydrocarbon components present in the fluid system (specifically how many "different" components are present) and whether δij has been adjusted extensively as part of the EOS adjustment.

In most cases, however, the ranking is low where many Eigenvalues are negligible or zero.

It is therefore assumed that an Eigenvalue in a size that is smaller than a certain drop tolerance of, for example |λα|< εtol, can be safely ignored without affecting predictions and where an effective ranking m ~ ranking (1 - δij occurs.

By approximate representation

can now be introduced and lead to significant savings in algorithms that are introduced later, provided that m«c applies.

It is worth emphasizing that 1- by is a constant matrix for a fixed fluid description and consequently the decomposition in (15) will only be calculated once and can be used in algorithms without incurring a runtime penalty

By replacing the approximate expansion (15) with (6), we get:

Then let M = m 1 and introduce the matrix Q(P,T) ∈ R<M>'<C>if the elements are defined by

With these definitions we can write (16) as

Introducing the vector (Q in ... QM) of reduced variables,

and substituting (19) for (18) gives a particularly simple form for A and its mole fraction derivative:

It also follows that

Note that the physical limitation is on mole fractions in the future, 0 < x ≤ 1 implies that the reduced variables are necessarily constrained by the min/max values of the matrix entries, i.e.:

It follows in particular from (23) that min(Bj)<B< max(Bj ) must apply.

Inserting equations (20)-(22) instead of A, ∑i and B in the original expression for the fugacity coefficient (14), the dependence is reduced from c components to the smaller M reduced variables,

where Z = Z(Q) followed by Z=Z(A,B) in light of equations (20) and (22).

3 Stability test algorithm in reduced variables

The stability of a mixture is determined by the tangent plane distance (TPD) criterion established by Michelsen, M.L, "The Isothermal Flash Problem, Part I. Stability", Fluid Phase Equilibria 9 (1982) 1-19 which can be briefly expressed as follows:

The phase of the composition z is stable at the given P, T and if and only if i

for any admissible sample composition y.

Global minimization problems of this type are not generally directly straightforward to deal with, but a reasonable compromise and performs a search for stationary points (25) and verifies non-negativity of the TPD at all such.

The stationary points of the TPD can all be found among the solutions of the equations

In yi+ In φi(y) Inzi -Ιηφi(z)= k (26)

and a necessary condition for stability is the requirement that k ≡ TPD(y)> 0 at all stationary points.

A suitable change of variables Yi ≡ yiexp(-k) transforms (26) into c equations in the unforced, non-normalized moles,

(27)

or, denote the particular part of (27) as dj≡ In z; In^j(z), rewritten more compactly

The necessary condition for stability is thus at all stationary points.

In order to carry out a comprehensive search for stationary points and ensure that no unstable states are missed, the standard procedure requires that the equations (28) be solved starting via both a "light" and a "heavy" test phase.

The first estimate in such calculations is based on Wilson K values provided that (T,P) by the expression:

InK<w>= In(Pci/P) 5.373 ( 1Wi)( 1 - Tc/T) (29)

3.1 Procedure for direct solution

This formulation, which through extensive testing has been preferred for a composite reservoir simulator, can be considered a Newton iteration applied directly to the system to define conditions for reduced variables in sample phase mols. This is carried out by using the provisions for stationarity (28).

To make the above concrete, one can consider the definition with reduced variables (19) written in residual form for the sample phase moles Y.

By introducing the condition of stationarity (28), and observing that the expression InΦi can be regarded as a function of Qa, by the expression the sample phase moles as a function of the reduced variables, i.e.:

(31)

Mixing the equations (30)-(3 1) into a closed system that must be satisfied by the M unknown, reduced variables Qaved a stationary point of the TPD function:

The Jacobin of the above system can be shown with the form

The derivatives of the fugacity coefficient in (33) can be calculated via EOS via (24)

After solving the Newton system,

the solution is updated via

(35) and the denormalized mols are then revised to meet TPD stationarity

reality (31), i.e.,

If the update (35) violates the bounds (23) on the RVs, the update is abandoned and a regular subsequent replacement step is performed instead. This is performed using step (36) with the old iteration Q<k>.

4 Phase split algorithm in reduced variable

When the proximity of two hydrocarbon phases has been determined (through e.g. stability analysis), the compositions and sizes of the equilibrium phases must be calculated, traditionally by solving sets of c non-linear equations with equal fugacity for e.g. the moles in the vapor phase.

In this section, the method with reduced variable is described as instead solving a set of (M 1) primary variables where P denotes a primary phase and β<p>is a corresponding phase fraction.

In principle, either phase can be named as primary and solved for, but from a numerical standpoint it is better to solve for the least abundant phase, especially near the phase boundaries. Thus will be used near a bubble point and near a dew point.

The least abundant phase is determined by solving the Rachford-Rice equation for vapor fraction based on the current feed composition ztog a prior guess for balanced K values, Ki ≡ yi/xi

The variables corresponding to the non-solution phase S (secondary) are calculated from a mass balance:

Note that it is not possible to divide by zero, assuming that the phase S is more correct.

Using a procedure outlined in Firoozabadi, A. and Pan, H., ''Fast and Robust Algorithm for Compositional Modeling: Part II - Two-Phase Flash", SPE 71603 solves the set of (M l ) equations (provided that steam is the primary phase),

Using the standard Rachford-Rice relationship to express this phase composition in K values gives

The system (39) is closed since the K values can be shown to depend only on primary variables. Note also that the Rachford-Rice equation (no. 2 in the set (39)) is not solved separately in this formulation, but rather part of the system of residuals.

4.1 Generalized K-value form

The system (39) must be linearized for X<p>for the least abundant phase. To continue, it will now be shown how the defined equations can be written generically for the P and S phases and "generalized" K values to make the phase switching practical.

The Rachford-Rice expressions using K values in regular form are:

By introducing reciprocal K values Ki= 1 /Ki it will appear that

This example makes it clear that by introducing generalized K-values that coincide with standard K-values when the primary phase is steam and otherwise equal reciprocal K-values, i.e.,

(40)

the secondary and primary phase compositions are given by

The phase splitting conditions can now be written uniformly as:

4.2 Jacobin system

Name the primary and secondary variables

Given the derivatives of the primary/secondary phase compositions with respect to the primary variables, the Jacobin system is easily assembled from the remaining definitions:

4.2.1 Derivatives of phase compositions

Using the formula for phase compositions in conjunction with generalized K values,

Insert, for example, K”'Ki -1 and we sign:

And

4.2.2 Derivatives of generalized K values

From the definition of generalized K values, it follows that

Using mass balance, primary variable derivatives of the secondary phase fugacity can be written as

4.2.3 Derivatives of fugasities

The derivatives of fugacity coefficients in connection with the RV of the corresponding

phase shown in equations (46) and (47) follows directly via (24).

4.3 Newton Update

After the construction of residuals (42) and composition of

The Jacobin system via equations (43)-(47), can the Newton system

is solved and new iterations are generated from

As previously mentioned in connection with stability testing, the updated, reduced variables must always fulfill constraint one (23). In this case, the phase shift must also be ensured to ensure that the physical limits 0 < β < 1 are respected. See fig. 3 which schematically shows these limits.

If violation of these constraints is detected at any point in the repeated sequence defined by (48), either the conditions given corresponding to a single phase or a poor quality initial guess rule out the convergence. In such cases, it is preferable to terminate the iterations and continue with flash calculations from the beginning.

5 Reservoir flash - A summary

Fig. 4 is a functional diagram of a non-linear repetition loop with flash calculations made during computer simulation of a fluid flow according to the invention in an underground, hydrocarbon-bearing reservoir model. Fig. 2 is a flow chart showing the combined use of the condition stability test and reduced variable transform to solve the flash problem at a particular time step of a composite reservoir simulator.

The PVT module in a reservoir simulator is always responsible for detecting the occurrence of new phases in each calculated grid cell since this information cannot be retrieved in any way to set more primary variables and reservoir equations.

For cells where there are equilibrium phases, the situation is different in that the corresponding equilibrium limitation appears in general systems of master equations,

It is possible to simply include the residual values (49) in the rest of the system, as this implies that simulator and flash residuals coincide without special flash calculations. The second possibility, known as "exact flash" determines equilibrium phases and compositions that meet (49) a strict tolerance. In the exact flash method, the residual values (49) become (numerically) zero due to performed flash repetitions. Other derivatives of (49) in conjunction with simulator variables are needed to take into account the contribution of equilibrium to the Jacobin simulator.

The exact flash method is preferred, partly because it leads to a modular structure, but mainly because the coincidence of highly non-linear flash constraints often requires special treatment. In addition, the extra work associated with performing accurate flash calculations will again be offset by a faster convergence of the simulator's non-linear iteration arising from more accurate equilibrium phases.

An effective application of stability testing of fixed-trial algorithms presented in Sections 3.1 and 4 to the iterative solution of flash problems encountered during the time steps in the simulator requires accurate use of existing conditions of pressure, temperature and composition on the simulation grid.

5.1 Two-phase area

In general, saturated cells are very likely to remain biphasic rather than transition to monophasic in subsequent iterations. In the vast majority of cases, equilibrium phase compositions and sizes that have already been calculated provide an excellent starting point for new phase split calculations at the updated cell conditions. Consequently, an attempt will always be made to proceed directly to the Newton procedure described in section 4.3 using previous K values as an initial guess and to converge the master equations (42) to a strict residual tolerance (default is Σsplit= 10<10> ) in a very small number of repetitions. Excessive iterations, or an excess in repeated variables, should this occur, lead to a quick termination of the iteration in favor of the first flash calculations.

5.2 Single-phase area

Cells that are undersaturated are more likely to remain single phase than to develop additional phase, but this can only be determined with absolute certainty through stability testing and the traditional method requires this testing to begin from a rough starting point, for example Wilson's K values, equation (29 ) and using a 'light' and a 'heavy' trial phase for each new condition encountered.

An appropriate purpose is to reduce the number of expensive stability tests that are carried out. This has been attempted in the past by limiting the testing to specific candidate cells, for example those bordering clusters of existing two-phase cells. However, such algorithms eventually become heuristic and introduce an iterative component in the sense that new neighbors must be tried after instability has been revealed. This has undesirable effects in the parallel simulator.

The preferred method associated with Claus P. Rasmussen, Kristian Krejbjerg, Michael L. Michelsen and Kersti E. Bjurstrom, Increasing the Computational Speed of Flash Calculations with Applications for Compositional, Transient Simulations, Society of Petroleum Engineers, SPE 84181, February 2006 SPE Reservoir Evaluation & Engineering and will hereby refer to Conditional Stability Testing (CST). Instead of considering local lattice conditions, this method is based on monitoring each individual phased cell's approximate location in the phase plane. The image below shows the cases that occur.

The single-phase area is further divided into a "shadow" area (C) next to the phase boundary and is characterized as the set of P and T for which TPD equations (26) have a non-trivial solution with a corresponding positive value of TPD. Outside the shaded region is the far region (D) where only a trivial solution (y = z) for the TPD equations exists. The shadow area thus acts as a buffer between the statuses far into the single-phase area and the two-phase area itself.

The central idea of CST is to try and skip the stability calculations in zone (D) provided that the magnitude of the change experienced in a given cell is sufficiently small. For the cell area (C), the stability test is "single-sided" and begins with the Newton iteration from the previous, calculated and non-trivial solution with positive TPD.

As the width of the shadow zone can be shown to shrink to zero in the critical region, a measure of the distance to the critical point becomes necessary to safely skip calculations in zone (D). This measure, introduced by Michelsen, is given by the smallest Eigenvalue of the matrix

When the conditions Ρ*, Τ*, zi* in a cell are found to lie in the area (D), the matrix (50) is calculated and its smallest eigenvalue λΒ= min(eig(B)) is determined. All values marked (*) are then saved for reuse.

In subsequent calculations under new conditions, P, T, z define; state variable changes relative to the base point which

(51)

A computational tolerance is determined by scaling a user-supplied tolerance parameter (default value εΑ= 0.1) by the tabulated, approximate distance from the critical point

The stability test can be omitted under new conditions, provided that the following conditions are all met:

5.3 Initial flash calculations

In the procedures outlined in sections 5.1 and 5.2 above, there are opportunities for algorithmic errors. For example, an assumed two-phase condition may actually be a single phase or the quality of the first guess insufficient to allow correction in Newton's method. Analogously, a single-phase state that was previously in the shadow or remote area may have changed in conditions that are too great to be able to draw the conclusion that the fluid remains stable at a single phase. In such cases, situations arise where there is no initial information, an initial calculation becomes necessary. The present procedure is based on Michelsen's work. The main difference lies in the use of the reduced, variable technique section 3.1 for the stability test step.

It should first be noted that for two phases of assumed composition x, y, the difference in GFE with the fluid can be considered as a single phase, expressed as:

By using mass balance in; = (1 - β)χί+ βyi; this can be described as

or in recognition of the tangent plane distance of liquid and vapor compositions,

ΔG = (1 - β)TPD(x ) βTPD{y)

If the compositions x, y can be found to satisfy the mass balance with AG < 0, the mixture becomes z; unstable at the current conditions P and T. If, in addition, either TPD(x) < 0 or TPD(y) < 0, the mixture also becomes unstable (in practice, a small threshold value is used instead of zero, the default setting becomes TPD z εstab< =>-10<11>). In such situations, the stability test becomes redundant.

The steps in the Ab Initio algorithm are as follows:

First evaluate φ ( P , T, z ), Gibbs Free Energy of the feed and Wilson K values

using equation (29). Then cycles of subsequent replacement are performed. If AG < 0.0 during this procedure is observed, the feed becomes unstable and current compositions can be used in subsequent split calculations. If instead it is the case that TPD(x) < 0 or TPD(y) < 0, the same conclusion applies. The choice of estimates will now be based on which phase indicated instability, i.e. if TPD(x) < 0, estimate K values as log If after three repetitions instability has not been revealed, no definite conclusion is possible and a full stability test must is performed.

If instability is detected, either through stability testing or the SS procedure outlined above, estimates now become available for which the objective function AG is negative. Three additional cycles each consisting of three SS repeats are used each followed by an attempt to accelerate the process. Only convergence to a minimum of GFE is possible with this method. If the convergence tolerance is met after completing the cycles, the logarithm terminates, if not, a second-order rigorous GFE minimization algorithm is used with strict descent applied for the final convergence. Consequently, the trivial solution is always avoided and convergence to a minimum of GFE is guaranteed.

5.4 Phase marking

After it has been confirmed that the fluid of composition z; is stable as a single phase at the current draw and temperature, a label of either 'oil' or 'gas' must be assigned to it, the main reason being the need to use the correct relative permeability table calculation of the flow properties of the phase.

However, it is important to highlight that such marking gradually becomes meaningless as the critical point approaches and that there is currently no universally accepted method for marking, as current simulators show a large variation in this regard. From a physical point of view, the flow properties of a phase cannot depend on this label as the phases cannot be separated from each other.

Although it will be seen that the most rigorous method of marking is the determination of the true critical point of the mixture, which is an expensive process and an accurate determination may not be necessary, especially if the properties are extra polished between oil and gas close to the true or approximate , critical point.

For these reasons, the present simulator uses a simple correlation for quasi-critical temperature which can be expected to be sufficiently accurate to correctly label the phase well away from miscible conditions. This so-called Li correlation represents a weighted average of the component's critical temperature,

Here, Γ is a correction factor typically unity unless the model needs to be adjusted to fit the initial data, i.e. the location of a gas-oil contact.

Using (52) at operating temperature T, a fluid that fulfills

labeled oil, or otherwise gas.

Fig. 4 shows a procedure for reservoir simulation. In step 100, the reservoir model and reservoir data are used. The fluid phase properties and equations of state are calculated in step 110. Such calculations can be as described in sections 1-4 here. A Jacobin matrix is generated in step 120 and may be as described in sections 1-4.2 herein. The linear equations are thus solved in step 130 according to the description in sections 1-4 — 2.3. The solution is then updated in step 140 according to the description in section 1-4.3. The solution is then driven for stability or convergence according to the description in sections 1-4.2. 3. Then the calculated solution is sent to the user in step 160.

In fig. 5, a method 200 for reservoir simulation is shown. In step 210, a cell is selected with a vapor phase and a liquid phase in the cell. A calculation is made whether the vapor phase and the liquid phase are present in a less abundant amount in step 220. In step 230, the phase with the least abundant amount is assigned as the primary phase and the other phase is assigned to the secondary phase. In step 240, the phase properties of the primary phase are calculated using the primary variables with the primary phase and the phase properties of the secondary phase are also calculated using the mass balance and the other variables associated with the secondary phase. In step 250, stability is ensured in the calculations by dividing by the primary phase rather than by a value close to zero.

In the method, the phase properties of the primary and secondary phases may include pressure, temperature and pressure of the primary phase. In the method, the phase properties of the primary and secondary phase may include pressure, temperature, pressure, composition and a size of the primary phase.

In the method, the calculations in step 240 may further include the following steps to calculate the phase characteristics of the primary and secondary phases: (i) using a reduced variable algorithm with the primary and secondary variables associated with the primary and secondary phases producing a Rachford-Rice expression, ( ii) linearize the Rachford-Rice expression with K-values and reciprocal K-values to thereby produce linear expressions, (iii) generate a Jacobin matrix by utilizing the primary and secondary variables and (iv) solve the linear expressions and in the Jacobin matrix to update properties and test for stability. Such steps are described in sections 1-4.2.3 here. In the above method, step (iii) may also include the steps: (1) calculate derivative phase compositions, (2) calculate derivatives of K values and (3) calculate derivative fugacity coefficients corresponding to each of the primary and secondary variables.

In fig. 2 and 6, another method 300 is shown for determining composition of fluid in a cell during a computer-activated reservoir simulation. In step 310, direct reduced variable split calculations are performed using K values when the cell had a multiphase fluid in a previous time step. In step 320, a single sided reduced variable stability test is performed using vapor incipient moles when the cell had a single phase fluid in the shadow region, i.e., the liquid side in the previous time step. In step 330, a single sided reduced variable stability test is performed using liquid starting moles when the cell has a single phase fluid located in the shadow region, i.e. the vapor side in the previous time step. In step 340, a first flash calculation is performed on the cell to determine fluid composition and on to step 360 when the cell is in the remote region. The first calculation can be as described in section 5.3 here. In step 350, it is determined whether there is an error in step 310, 320, or 330. In step 350, a first space calculation is also performed on the cell to determine fluid composition and then the process continues to step 360 when it is determined that an error exists. In step 350, when the fluid is determined to be a single phase, additional calculations are performed to determine the location in the phase plane and steps 320-340 are repeated. In step 360, the calculated results are used when there are no errors.

In fig. 7 in the system 390, a computer readable medium 410 is shown used during a reservoir simulation to determine composition of the fluid in a cell. As will be apparent to one skilled in the art, the computer-readable media may also be a component of a system where they form machine-readable media or software 410 to interface with an input device 400, such as a computer terminal and a central processing unit (CPU) 460. A person skilled in the art will also easily understand that such a system can also be part of a computer network. The computer media 410 includes a data receiver 420 that receives reservoir model signals and data from a source. The computer media 410 also includes at least one abundance allocator 430 that calculates whether a vapor phase and a liquid phase of fluid in the cell are present in at least an amount in response to input data received by the data receiver. The excess amount allocator 430 also assigns the phase having the calculated least abundant amount as the primary phase and assigns the other phase as the secondary phase.

The computer media 410 also includes a fluid phase property calculator 440 which calculates phase properties of the primary phase by utilizing primary variables in the primary phase. The fluid phase property calculator 440 also calculates phase properties of the secondary phase by utilizing mass balance and secondary variables associated with the secondary phase. The fluid phase property calculator 440 thereby ensures stability by dividing by the primary phase rather than dividing by a value close to zero. The computer media 440 also includes an output producer 450 adapted to produce and communicate the calculated phase properties of the fluid to a readable format, for example to a screen or a printer.

In the computer readable media 410, the fluid phase property calculator 440 may also include a reduced variable algorithm module 470 that produces a Rachford-Rice expression with the primary and secondary variables associated with one of the primary and secondary phases. The fluid property calculator 440 may also include a linearization module 480 that produces linear expressions from the Rachford-Rice expression with K values and reciprocal K values. The fluid property calculator 440 may also include a Jacobin matrix generator 490 that generates a Jacobin matrix by utilizing primary and secondary variables. The fluid property calculator 440 may also include a solution and stability tester 500 that solves the linear expressions and the Jacobian matrix to update the phase properties and tests for stability. In the Jacobin matrix generator 490, a phase composition derivative submodule that calculates the derivative of the phase compositions, a K value submodule that calculates derivatives of K values and a fugacity submodule that calculates derivatives of fugacity coefficients corresponding to each of the primary and secondary variables can also be used in the Jacobin matrix generator 490.

In fig. 2 and 6 another method 500 is shown for determining composition of fluid in a cell during a computer activated reservoir simulation. Step 510 determines whether the cell had a single phase or multiple phases in the previous time slice. From step 510, other steps are performed depending on the determinations in step 510. Step 520 is performed if the cell had multiple phases in the previous time slice. In step 520, direct reduced variable split calculations are performed using K values. Step 530 is performed if the cell had a single phase in the previous time step and was in the shadow region, ie the liquid side of a phase plane. In step 530, a single-sided, reduced wire at the stabilization test is performed using vapor-initiating moles. Step 540 is performed if the cell had a single phase in the previous time step and was in the shadow region, i.e. the vapor side of the phase plane. In step 540, the single sided reduced variable stability test is performed using liquid starting moles. Step 560 is performed if the cell is in the far region of the phase plane. In step 560, a first calculation on the cell is performed to determine fluid composition. In step 560, and after performing the first calculation, the method continues to step 580.

Step 570 is performed after performing steps 520-540 to determine if there is an error in fast processing. If there is an error, then an initial flash calculation is performed on the cell to determine the fluid composition. If the fluid is found to be a single phase, additional calculations are performed to determine location in the phase plane for use in subsequent iterations. Step 580 is then performed if there are no errors. In step 580, the calculated results are used.

Although the invention has been described above in connection with certain preferred embodiments and many details have been presented for illustration purposes, it will be clear to a person skilled in the art that the invention can be changed and that certain other details described here can vary significantly without deviating from the main principle of the invention.

For example, the invention may also include a system and a computer-readable medium that executes instructions for performing a complex reservoir simulation of an underground hydrocarbon reservoir. This system of computer hardware and storage will be able to perform the reservoir repacking procedure outlined above. Likewise, computer-readable media can have instructions for carrying out a complex simulation of an underground shy dr ok arb onreservoir in accordance with the aforementioned principles.

Claims (12)

P a t e n t k r a v 1. Procedure for complex reservoir simulation of a fluid flow, characterized by: (a) providing a reservoir model for an underground reservoir wherein the model has a grid defining a plurality of cells in such a manner that each cell is associated with the fluid therein, the fluid being in a vapor phase and a liquid phase, (b) selecting at least one cell of the plurality of cells associated with fluid that is in the vapor phase and fluid that is in the liquid phase (c) estimating which of the vapor phase or the liquid phase is present in the least abundant amount; (d) allocating the least abundant phase as the primary phase and allocating the other phase as the secondary phase; (e) calculating phase properties of the primary phase responsive to primary variables associated with the primary phase and calculating the phase properties of the secondary phase using the mass balance and responsive to secondary variables associated with the secondary phase, so that stability is ensured during calculation of the phase properties of the secondary phase due to allocation of the phase with the least abundant amount as the primary phase, (f) outputting a visual display responsive to the calculated phase characteristics of the primary phase and to the calculated phase characteristics of the secondary phase. 2. Method according to claim 1, characterized in that the phase properties of the primary phase include pressure, temperature and composition and quantity of the primary phase. 3. Method according to claim 1, characterized in that the phase properties of the secondary phase include pressure, temperature and composition and quantity of the secondary phase. 4. Method according to claim 1, characterized in that step (e) further comprises the following steps to calculate the phase properties of the primary and secondary phase: (i) use a reduced variable algorithm with the primary and secondary variables associated with the primary and secondary phases producing a Rachford-Rice expression, (ii) linearize the Rachford-Rice expression with K-values and reciprocal K-values to thereby produce linear expressions , (iii) generate a Jacobian matrix using primary and secondary variables, and (iv) solve the linear expressions and the Jacobin matrix to update the phase properties and test for stability. 5. Method according to claim 4, characterized in that the step (riding) further comprises: (1) calculate derivative phase compositions, (2) calculate derivatives of K values, and (3) calculate derivatives of fugacity coefficients corresponding to each of the primary and secondary variables. 6. Procedure for determining composition of fluid in a cell during computer-activated reservoir simulation, characterized by the steps: a) providing a reservoir model for an underground reservoir where the model has a grid defining a plurality of cells in such a way that each cell is associated with the fluid therein; (b) performing one of the following calculations for each cell based on the composition of the fluid in the cell at the previous time step; (i) performing direct reduced variable phase split calculations using K values when the composition of the cell has a multiphase fluid in a previous time step, where the direct reduced variable phase split calculations are performed such that the phase composition of the fluid present in at least abundant amount is the primary phase (ii) performing a single-sided, reduced, variable stability test using vapor-initiating modules when the composition of the fluid in the cell comprises a single phase placed in the shadow region, i.e. the liquid side in the previous time step, (iii) performing a single-sided, reduced, variable stability test using liquid-initiating modules when the composition of the fluid in the cell comprises a single phase located in the shadow region, i.e. the vapor side in the previous time step, (iv) not to perform a stability test when the composition of the fluid in the cell comprises a single phase in the phase region at the previous time step, (c) to determine whether there is an error in step (b) and to perform one of the following steps: (i) performing a first flash calculation to determine fluid composition in c lien when it is determined that there is an error in step (b); (ii) using the calculated result from step (b) when determined not to be in error in step (b), and (d) outputting a display representing the composition of the fluid in response to step (c). 7. System that is braked during reservoir simulation to determine the composition of fluid in a cell, the system being characterized by including: a data receiver that receives input data from a source, where the input data comprises a reservoir model and data, a processor, a computer program that can be executed on the processor, where the program includes: (a) a least abundant amount of allocation module which (i) calculates whether a vapor phase or a liquid phase of fluid to the cell is present in at least abundant amounts in response to input data received by the data receiver; and (ii) assigns the phase that has the calculated least abundant | as the primary phase and assigns the second phase as the secondary phase, (b) a fluid phase properties calculator which (i) calculates the phase characteristics of the primary phase responsive to the primary variables associated with the primary phase, and (ii) calculates the phase properties of the secondary phase using mass balance and secondary variables associated with the secondary phase, the fluid phase properties calculator ensuring stability during calculation of the phase properties of the secondary phase due to the least abundant allocation module assigns the least abundant phase as the primary phase, and (c) an output producer adapted to produce and communicate the calculated phase properties of the fluid in a readable format. 8. System according to claim 7, characterized in that the fluid phase property calculator includes: a reduced variable algorithm module that produces a Rachford-Rice expression with the primary and secondary variables associated with the primary and secondary phases, a linearization module that produces linear expressions from Rachford-Rice expressions with K values and reciprocal K values, a Jacobin matrix generator that generates a Jacobin matrix responsive to primary and secondary variables, and a solver and a stability test that solves the linear expressions and the Jacobin matrix to update the phase properties of the primary and secondary phases and tests for stability. 9. System according to claim 8, characterized in that the Jacobin matrix generator further comprises: a phase composition derivative submodule that calculates derivatives of primary and secondary phase compositions, a K-value submodule that calculates derivatives of K-values, and a fugacity submodule that calculates derivatives of fugacity coefficients corresponding to each of the primary and secondary variables. 10. Method for determining composition of fluid in a cell during a computer-activated reservoir simulation, the method comprising the steps: (a) providing a reservoir model for an underground reservoir wherein the model has a grid defining a plurality of cells such that each cell is associated with the fluid therein; (b) determining whether each cell has a single phase or multiple phases in a previous time step; (i) if the cell had multiple phases in the previous step, perform directly reduced variable phase split calculations by selecting K values, (ii) if the cell had a single phase in the previous step, perform the following steps depending on whether the cell's previous position in a phase plane: (1) if the fluid is in the shadow region, i.e. the liquid side, perform a single-sided reduced variable stability test using vapor-initiating moles, (2) if the fluid is in the shadow region, i.e., the vapor side, perform a single-sided, reduced variable stability test using liquid-initiating moles, (3) if the fluid is in the remote area, do not perform a stability test of the cell to determine fluid composition and then skip to step (c), (c) determine if there is an error by rapid processing at (i), (ii) or (iii) ), (i) if there is an error, perform an initial flash calculation on the cell to determine fluid composition, (1) if the fluid is found to be a single phase, perform additional calculations to determine location in the phase plane for use in subsequent iterations; (2) if the fluid is found to be in a plurality of phases proceed to step (d) (ii) if there are no errors, proceed to step (d); (d) emitting a display representing the composition of the fluid. 11. Computer program stored on a computer-readable medium used for complex reservoir simulation of fluid flow in an underground reservoir, where the program includes: a least abundant allocation module which estimates which of the fluid in the vapor phase and the fluid in the liquid phase is present in the least abundant amount in a lattice cell, and assigns the phase with the least abundant amount to the primary phase and assigns the other phase to the secondary phase, a fluid phase properties calculator that calculates phase properties of the primary phase that is responsive to the primary variables associated with the primary phase and calculates phase properties of the secondary phase using mass balance and that is responsive to secondary variables associated with the secondary phase, and an output producer module that communicates the calculated phase properties in a readable format to a visual display. 12. Computer program according to 11, wherein stability is ensured while the fluid phase properties calculator calculates the phase properties of the secondary phase because the least abundant allocation module assigns the least abundant phase as the primary phase.