EP2030147A2 - Efficient application of reduced variable transformation and conditional stability testing in reservoir simulation flash calculations - Google Patents

Abstract

Methods and computer readable media for performing a compositional reservoir simulation of a subterranean hydrocarbon-bearing reservoir is provided (500). Reduced variables for flash computations are utilized with a methodology of conditional stability testing for achieving optimal efficiency of phase behavior computations in a compositional reservoir simulator (510). A least abundant phase is selected as primary variables for a primary phase and a secondary phase is selected for a more abundant phase to ensure stability by not dividing by a value near zero due to the selection of the primary phase as being associated with phase which is the least abundant. A bounded interval may be used to limit solution changes in reduced variable algorithms to achieve great stability of algorithms. Stability test may be performed During flash computation using reduced variables by employing Direct residual from based on the definition of the reduced variables and the tangent-plane distance condition.

Description

EFFICIENT AJ¹PLICATION OF REDUCED VARIABLE

TRANSFORMATION AND CONDITIONAL STABILITY TESTING IN

RKSKRVOIK SIMULATION FLASH CALCULATIONS

Ii^IIEICALfIELD

The present invention relates generally to computer enabled reservoir simulation of fluid flow in subterranean reservoirs,, and more particularly, to compositional reservoi r s irn ulation.

BACKGROUND Of. ΪHE MYMXfiibl

Conceptually a compositional reservoir simulator for simulating How in a subterranean hydrocarbon-bearing reservoir can be viewed as modeling a series of connected mixing tanks of fluids (cells) at given pressure, temperature and compositions. As time evolves (as the simulator is taking time-steps toward some final time at which results are sought), conditions in the tanks change as a result of fluid movement, wells and other external factors. Flash, calculations are necessary to establish, for each new set of pressure, temperature and overall fluid composition, the number of .fluid phases, their amounts and compositions. This calculation fundamentally involves finding the minimum of a thermodynamic state function (the Gibbs Free Energy (GFF), and as such is iterative in nature and frequently difficult to converge and computationally expensive, particularly when detailed Quid models are used, i.e., when there are many hydrocarbon components present. There is therefore great interest, in. devising algorithms which are computationally efficient yet robust and accurate.

for the purpose of discussion, the activity performed in "flash" calculation shall be subdivided into stability testing, which attempts to reveal instability of a given phase at the current conditions; and split calculations, which aims at determining the eςμiiHbrium phases and compositions for an. assumed phase configuration. One approach to increasing the efficiency of flash calculations In a computational reservoir simulator is described by Claus P, Rasmussen, Kristian Krejbjerg, Michael I... Miehrfsen and Kersli E. Bjuratrom, Increasing the Computational Speed of Flash Calculations with Applications for Compositional, ^'Transient Simulations, Society of Petroleum Engineers, SPE 84181, February 2006 SPF Reservoir Evaluation & Engineering, Tn performing flash calculations, the majority of time is spent doing stability analysis. Rasmiissen et al, proposed criterion for bypassing many of the stability analysis cheeks.

Referring to FlG.1. a pressure-temperature map is shown for a fluid in a cell. Point A is shown in a two-phase region where both a gas phase ami a liquid phase exist. Point B is located on a transition line between a two-phase region and a one phase region (the phase boundary). Point C is located in a "shadow zone" of the one-phase region, close to the two-phase region. Finally, point D lies in a ''remote'^* region far into the single-phase domain, Also, a vertical line is shown which separates single-phase liquid on the left and on the right is single-phase gas. Depending on where the estimated phase state of fluid in a cell, certain stability calculations may be omitted rather than performing stability analysis for ail cells during iterations in a time. This general criterion for bypassing calculations shall be described in greater detail below with respect to the stability testing in Section 5,

The sub-steps corresponding to "split" and "stability^" calculation described in Rasmussen et. al. use a "iraditiona!" approach, with the attendant solution of nonlinear problems of size equal to the number of hydrocarbon components. There is scope for improving the efficiency of these sub-steps, particularly for simulation models involving a large number of components.

Firoozabadl, A. and Pan. 11.. Fast and Robust Algorithm for Compositional Modeling; Part ϊ Stability Analysis, SPE 63083 and Firoozabadi, A. and Pan. H., Fasi and Robust Algorithm for Compositional Modeling: Pan 11 ~ ^'I'wo-Phase Plash, SPR 7160? discuss the application of reduced variable strategies for stability and split calculations in compositional reservoir simulation: however the authors do not teach how stability tests can be avoided. Also, the particular stability algorithm, as formulated, may experience convergence difficulties, particularly when encountering conditions far into the undersaUirated zone, In addition, the split algorithm is formulated in terms of the vapor phase and will exhibit numerical and/or convergence difficulties near dew-points, due to the virtually non-existent liquid phase.

Newton^'s method is commonly used in solving nonlmeai systems of equations. Care must be taken to ensure that iterates do not exceed physical bounds on the unknowns. In applying Newton^'s method to problems in phase behavior formulated in terms of reduced variables, there is a need to ensure that physical bounds on the reduced variables are not violated.

The shortcomings of previous methods for compositional reservoir simulations cited above will he addressed by the detailed description of the invention that follows below.

SUMMARY OF THE INVENTION

A method, system and computer readable media carrying instructions for performing a compositional reservoir simulation of a subterranean hydrocarbon- bearing reservoir is provided. Reduced variables tor Hash computations are utilized, combined with a methodology of conditional stability testing for the purpose of achieving upttmai efficiency of phase behavior computations in a compositional reservoir simulator.

Abo, preferably a least abundant phase is selected as primary variables associated with a primary phase and a secondary phase is selected for a more abundant phase such that stability is ensured by not dividing by a value near zero due to the selection of the primary phase as being associated with phase which is the least abundant. Also a bounded interval may be used to limit solution changes in reduced variable algorithms (phase split and stability) to achieve greater stability of algorithms. Further, stability testa may be performed during flash computations using reduced variables, by employing a direct residual form based on the definition of the reduced variables and the tangent-plane distance condition.

~ 1 - It Is an object of the present invention to combine the concept of reduced variables for flash compulations with a methodology of conditional stability testing for the purpose of achieving optima! efficiency of phase behavior computations in a compositional 5 reservoir simulator.

It is an object of the present invention to provide a. more reliable reduced-variable phase-spin algorithm by selecting primary variables eoπespoπding to the least abundtuH phase present; 10

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

i τ> It is yet another object to provide an enhanced method for performing stability tests using reduced variables by employing a direct residual form based on the definition of the reduced variables and the taogeni-piaπe distance condition,

BRIEF¹ _. DESCRIPTION OP I Ht- DRAWINGS 0

These and other objects, features and advantages of the present invention will become bεttei understood with regard to the following description, pending claims and accompanying drawings where:

25 Vi.il 1 is a pressure-temperature diagram delineating several regions in the phase- plane to illustrate concepts centra! to a conditional stability teat approach;

FϊG. 2 is a flow-chart illustrating the combined usage of conditional stability lest logic for the overall flash update, and reduced variable algorithms for iterative solution of 30 the stability and phase-split problems at a particular time-step of a compositional reservoir sim uiator: FIG. 3 illustrates physical limits applying to the reduced-variables for the purpose of safe-guarded Newton iterations:

FlG. 4 is a functional block-diagram of a non-linear iteration loop including flash 5 calculations made during computerized simulation of fluid flow, incorporating the present invention into the context of a subsurface hydrocarbon -bearing reservoir model ;

FIG. 5 is a functional block-diagram of an embodiment of a method in accordance 0 with the present invention;

.FlCi. 6 is a functional block-diagram of another embodiment of a method in accordance with the present invention; and

5 HG. 7 is a schematic representation of an embodiment of a system and computer readable media in accordance vAih the present invention,

DETAIlI-D DESCRIPTION OF THE IN VHNTjQN

0 The following nomenclature shall be used with the equations that follow:

Sχmbo|s ^β ~ Tangent Plane Distance (TPD) f -^{? ~} Gibbs Free Energy function (GFC) 5 ^{c :::} Number of hydrocarbon components ^m ------ Number of non-zero Eiger. values

M - Number of reduced parameters. Equals ^{m +' J}

^ ^::: Pressure

^ === Reduced variables (M-vector) of entries ~"

?0 Q = Reduction coefficients matrix of size ^jWc of entries^"^' i? ^~ Universal (las Constant T Temperature x L iquid phase corn position (mole fractions} y Vapor phase composition (mole fractions)

Y Unnormali/ed moles of trial phase Feed (lota!) composition (mole fractions)

Subscripts

Component index a Reduced variable index

L Liquid

V Vapor

Superscripts

F Feed phase

T Trial phase

Greek symbols

^>f>< = Fugacity coefficient

Binary interaction coefficients (BIC;

The teachings of Fϊroozabadi. A. and Pan, FL, Fast and Robust Algorithm for Compositional Modeling: Part / - - Stability Analysis, SPE 63083 and Firoozabadi, A. and Pan, H., Fast and Robust Algorithm for Compositional Modeling: Part II — Two- Phase Flash, SPE 7160, are hereby incorporated by reference in their entireties. Similarly, the contents of U.S. patent application, 2006/0036418 Jo Pita et a!.. Highly- Parallel, implicit. Compositional Reservoir Simulator for Multt-MiHion-Cell Models, is incorporated by reference in its entirety. Further, the contents of Clans P. Rasmussεπ, Kristiaπ Krejbjerg, Michael L, Micheisen and Kersti E. Bjurstrom, Increasing tha Computational Speed of Flash Calculations with Applications for Transient Simulations, Society of Petroleum Engineers, SPE 84181, February 2006 SPL Reservoir Evaluation & Engineering, are incorporated reference. Finally, the teachings contained within Michael L, Michelsen, The hoihtπnal flash problem, Part I. Stability, Fluid Phase Equilibria, 9 ( 1982) 1 -19 Michael L. Michelsen, The hot kmna! Flash Problem, Part IL Phase-split Calculation, Fluid Phase Equilibria, 9 (1982; 21 -40 are also incorporated by reference in their entireties.

HG. 4 show- the genera! steps taken during a non-linear iteration loop. Property and EOS calculations are made, A Jacobiau matrix is then generated. A linear solver is used to solve a linear set of equations for a solution. The solution is then tested for sufficient convergence. W not sufficiently converged then new EOS and property calculations aie made Otherwise, the converged results are output.

1 Cubic Equation of State

For a pure substance an Equations of State (EOS^') is a mathematical relationship between pressure, temperature and volume; for a mixture, composition is added Io thi^ relationship. The cubic form of the EOS is by far the most popular, and, in particular. {he Redlich-KwoΩϋ.~Soave-Peng-Robirtson family of EOS has Song been {he industry standard in compositional reservoir simulation.

The preferred BOS is written geneήeally in the pressure-explicit form

to encompass all members of this family.

In (1), parameters »*, andm, are used to represent the influence of the temperaiure- depeαdeπt attractive terra a ■^■-^■-^■ αC/'jand ^ is the repulsive term.

The specific FX)S corresponding to parameter selections rø, and m, is given below:

Defining the compressibility factor, or Z-facior.

and introducing non-dimensioπai counterparts A, B to EOS parameters α, έ- through expressions

υ , = — nπr (4) a re-arrangement of ( I) using the definitions (2}-f4) yields the cubic form υf the EOS,

.S which is solved depending on phase ami other considerations for the appropriate rooi Z = Z(A, B).

When applying the EC)S to a mixture - as opposed to a pure substance - mixing rules i 0 are appl ied to calculate the parameters a and b (or, equiva JentJy, A and 5 ).

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

Here, ι},. are the Binary Interaction Coefficients (^"BIC), accounting for chemical smeraeiions between components of dissimilar type. It is a symmetric matrix with zero diagonal entries, 20

The repulsive parameter is customarily calculated as a molar average,

25 ϊn the expressions (6 )-(?), the pure-component terms arc given by;

⁽ S) Ω, E i) - -^- (9)

which introduces reduced pressures and temperatures, P, ≡ P / P₁. and T_^ ≡ Trl] , respectively, where I]. and /' denote component critical properties.. The functional form of the temperature dependent function Ω.₍ , - Ω_:, . [T) depends o« the ROS chosen. Using >r. to designate component asceittric factor, the defining rcuttionships are as ^"follows:

, -i?

Ω (T) .: Q^ j 14 (0.48 ^•!^• 1 ,574w, -0.176 wf )( 1 -^T. i \ for SRK

Ω (T) - Q H<0.3746Φt l .54226w₁ -O.26992w? K l-,/7: Η lor PR (12)

The "corrected" form of Peng-Robinson is also supported, providing an alternative to (12) of Ae form

> 0.49

(13)

The term D_j , ~ Ω, is a constant mdepeπdεnt of teniperatυre for the Redlich-Kwong- Soave-Perig-Robiϊiaon family of EOS. The defaυli values of she EOS constants Ω._r and Ω_;, are given by the table below - they can be overridden by ihe user:

KOS Ω Ω,

RK.SRK 0.4274802 0.08664035 PR 0.45723S529 0.07796074

A traditional weakness of two-parameter EOS such as (! ) is frequently poor prediction of liquid density. To remedy this shortcoming, a three -parameter extension; via the standard Peneioux at. al, volume shifts may be used, Jn this framework, the molar volume is calculated according to

V " V

Here, r ^'"^'■" is the molar volume predicted by the EOS, equations (2) and (5); and th. correction term is calculated from

v- = V λ'.C,

where :%, detxMes phase composition and C₁ is a set of volume-shifts, related to user- supplied nun-dimensional volume shifts,?.;, according to

^■ sp, ^■-^■-^■ s. — y

Finally, iiigaeity coefficients and their derivatives are fundamental building blocks m the construction of EOS algorithms. They can be calculated directly irom ihe BOS using first principles. For an EOS of the generalised type U) they can be show to have the form

4 ?y. B. , ,Z + m,B B. _{( ry}

-In(Z - S) 4 ^' t-^-ⁱ 1) {14}

Im, - HI, ) B A

The fogacuy of a component is a measure of its '^'tendency to escape" from a phase and therefore directly useful in equilibrium calculations, Hie condition of thermodynamic equilibrium can be expressed as/'' - f; ^' , or equivalent!}', x,φ; - y,φ* , from which the definition of K-vaSues can be introduced as

2 Reduced Variable Approximation to the EOS

As noted above that the matrix l ~ d^" _;i is symmetric and consequently endowed with a full, orthogonal set {λ_{(i >}vf^' },a ~ \ ... c oϊ eigenvectors and coixespooding real Ci gen values. By the spectral expansion theorem this matrix can therefore be represented as a series

ⁱ -S^, = ∑V>?

The actual iank of the matrix 1 - <5_I; wiU depend on the number of nαn-hydrocarbon components present in the fluid system (specifically, how many '"dissimilar" components are present) and whether ilκS, have been adjusted extensively as part of the EOS tuning process.

However, hi the great majority of cases, the rank is low, with many Bgen values negligible or zero.

Assuming, therefore, that any Eigeu value in magnitude less than o certain drop tolerance, say |Λ,, | <&\,_{; !} can be safely ignored without affecting predictions an effective rank m * rankil - S₁₁ }: results,

An approximate representation,

(15) can now be introduced, which results in considerable savings in algorithms that will be introduced subsequently, provided thai m « c holds.

H is worth emphasizing thai 1 - ^_,, is a constant matrix for a fixed fluid description; hence the decomposition inherent in ( 15) will be calculated only once and can be used in algorithms without incurring a run-time penalty.

Substituting the approximate expansion (1 5} into {6) produces:

Next, IeU/ - ?κ + hind introduce the matrix Q(P, rj e PJ'^' ύm elements of which are defined bv

With these definitions, we can write (16) as

Introducing the vector (Q, , ... g_u ) of reduced variables,

α ^s ∑9»Λ ^{for α si 1}'- -^M ( 19) and substituting 09} into (J 8), yields a particularly simple form for A &nά its mole- fraction derivative:

ft also follows thai

For future reference, note thai the physical constraints on the mole fractions, O ≤ x_, ≤ l _> implies that the reduced variables are necessarily constrained by the miii'max values of the matrix entries, i.e.,

Ii follows in particular from (23) that rnin(i?, } < B < max( iS"₍) must hold.

Substituting equations (20)-(22) for /U ∑, and B into the original expression for fugacity coefficient (14), the dependency is reduced from c components to the smaller set oϊ^'M reduced variables.

wht-reZ = Z(O) follows (rom Z = 2(.4. m 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. !...„ 'The Isothermal Mash Problem, Part I, Stability^"', Fluid Phase Equilibria 9 ( 1 982) 1 -19 which can be stated succinctly as follows.:

The phase of composition - is stable at {he specified PJ^' if and only if

TPDi y) si ^T V₁ (In v, -t- In $ ( v) - In ;;, -- Inφ, (::}} > G C25)

/or any admissible (rial compvsidori y

Global minimization problems of this type are not generally directly tractable: however, a reasonable compromise is to perform a search for stationary points of (2S) and verity non-negativity of the ^"FPD at all such points.

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

In >', H- in φ_t i y) - in r , - In φ_t ( Ξ ) ~ A' (26)

and a necessary condition for stability is the requirement all stationary points.

A convenient change of variables V₁ s y, exp{-A) transforms (26) into c equations in the unconstrained, unnαπnaiized moles,

or, denoting the constant part of (27) as d. ≡ In 2. t- In φ_t (z) , re-written more compactly as

b .; + ln^ (J^'} -- 4 - 0 (28!

- 1 S ~ The necessary condition for stability is thus aϊ ail stationary points.

To make an exhaustive search for stationary points and ensure that no unstable states are missed, the standard procedure requires equations (28) to be solved starting from both a Might' and a "heavy' trial phase.

The initial estimate in such calculations is based on the Wilson K-values, given at ( T, P) by the expression :

3.1 The Direct Solution Approach

TbIs iorcnuiation, which through extensive testing has become the preferred for a compositional reservoir simulator, can be viewed as a Newton iteration applied directly to the system of defining relations .for the seduced variables in terms of rπal- phase moles. This is done while using the conditions of itationarity (28).

To make the above more concrete, consider the definition of reduced variables, i^'19), written in residual form for \hs tnal-phase moles Y ,

Invoking the condition of statkmarity, (28), and observing that the term In ^, can be viewed as a function of ^ , we can express the trial phase moles as functions of the reduced variables, i.e.,

Taken together, equations (3O)~{3 ! ) form a closed system that must be satisfied by the Λ/ unknown reduced variables Q_n &l any stationary point of the TPD function:

The Jacobian of the above system can be shown to have the form

The derivatives of the fugacity coefficient appearing in (33) can he caicυlated from the EOS via (24 )

VQ\ lowing solution of the Newton system,

the solution is updated via

and the tiiinoπnaiized iπoies are subsequently revised to satisfy the TPD sUitkmarity relations (31), i.e..

Lf the update (35_.) should violate the bounds (23 ) on the RVs, the update is abandoned and a conventional successive-substitution step is peifoπrted in its stead. This is accomplished by applying the step (36) with the old iterate Q^"' .

„ ^■ 7 „ 4 Phase Split Algorithm in Reduced Variables

When the existence of two hydrocarbon phases has been established (through e.g. stability analysis) the compositions and amounts of the equilibrium phases must be 5 calculated, traditionally by yoking a set of a nonlinear equaf-fugadty equations ϊov e.g. the moles in the vapor phase.

In this section, a reduced variable approach is described which allows us to instead solve a set of {M t- 1) primary variables X^1' -~ {Q{\...J2{: - β '^' ) ■ where .P designates a K) Primary phase and β^r is the corresponding phase fraction.

In principle, either phase could be designated as primary and solved for: however, from a numerical standpoint if is better to solve for the hast abundant phase, in particular near phase boundaries;. Thus (Q, ,... Ql, , V) would be used near a bubble-

! 5 point and (O₁'^' , .. , ⁽j{) , L ) near a devv-pυint.

The least-abundant phase may be determined bv solving the Rachfovd-Riee equation for vapor fraction, based on the current feed composition z, and a previous guess for equilibrium K-values. K. ≡ y. / x, .

20

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

15

Note that there is no possibility of division by zero, provided the assumption that the phase S is more abundant is correct. Uιili/.ing a methodology outlined in Firøozabadi. Λ. and Pan. I L "^"Fast and Robust Algorithm for Compositional Modeling: Part 11 - Two-Phase Flash", SPH ? 16(33 solves the set of (Λ/ -t- Ij equations (assuming Vapor is the primary phase),

Using standard Rachford-Rice relations to express phase compositions in terms of K- values yields

The system (39) ia closed, since K-values can be sho^"wn to depend only on the primary variables. Note also that the Rachford-Rice equation (2^ao in the set (39)) is fi'of solved separately in this fbmiυfation, but rather becomes pan of the system of residuals.

4. i Generalized K-value Form

The system (39) must be linearized in terms, oLY'^' for the kast-abuodaut phase. To proceed, we now illustrate how the defining equations can be written gencrically in terms: of i* and ά' phases and "generalized'' K-values. to make the phase-switching more convenient.

The Rachford-Rice expressions using Revalues in the customary form are:

Introducing reciprocal K-values K, s? 1 / /C, we see that

^'This example makes clear that by introducing generalized K-valυcs which coincide with standard K -values when the primary phase is Vapor, and otherwise equal reciprocal K- values, i.e.,

the secondary- and primary phase compositions are given by

! S where

The phase-split conditions can now be written uniformly as: 0

4,2 System Jacobian

Designate the primary- and secondary variables A" -- (££ , Q_\'-_{: >}β' _> ■

Given the derivatives of primary ^/secondary phase composition's with respect to the primary variables, the system jacobian is easily assembled from the residual definitions:

4.2.1 Derivatives of Phase Compositions

Recall formulae for the phase compositions in terms of generalized K- values,

Setting for convenience K; K₁ - 1 and v, -- , we find

⁽ξ:Υ

and

4.2.2 Derivatives of Generalized K-values

From the definition of generalised. K-vaJaes. K₁ ≡ φ^^'1 fφ." U follows thai

Using mass-balance, the primary variable derivatives of the secondary-phase iυgacUy tiffl be written

4.2.3 Derivatives of F υgacities

The derivatives of lligacity coefficients with respect to the RVa of (he corresponding phase appearing in equations (46) and (47) follow directly from (24).

4,3 Newton Update

Following construction of the residuals (42) and assembly of the system Jacobiaπ through equations (43}-(47), {he Newton system or &x^p-^k (48 ) lar

can be solved, and new iterates generated from

As mentioned previously in conjunct ion with stability testing, the updated reduced variables must always satisfy the constraints (23) . In the present case, the phase fraction JC^₊, « /T must also be safeguarded tυ ensure that the physical bounds 0 < /3^{' J'} ≤ l are respected. Sec FiG. 3 which schematically shows these bounds.

If violation of these constraints should be detected at any point in the iterative sequence defined by (48), either the conditions specified correspond to single-phase, or a poor-quality initial guess precluded convergence. In such eases, it is preferable to terminate the iteration and continue with ab inito Hash calculations.

5 Reservoir Flash - A Summary

FTG. 4 is a functional block-diagram of a non-linear iteration loop including Hash calculations made during computerized simulation of fluid flow according to the present invention in a subsurface hydrocarbon- bearing reservoir mυιk'1. FlG. 2 is a Dow-chart illustrating the combined usage of conditional stability test and reduced variable transformation for solving ihc flash problem at a particular time-step of a compositional reservoir simulator.

The FVT module in a reservoir simulator is always responsible for detecting the emergence of new phases in each computational grid ceil, since this information cannot in any way be inferred from the set of primary variables and reservoir equations For ceils in which coexisting equilibrium phases exist, the situation is different in the sense that a corresponding equilibrium constraint appears in the overall system of governing equations,

5 r^l * x,$^{' (}PS.x) ^~ y$^'i PJ\ y) for i ~ l... < c (49)

it is possible to simply incorporate the residuals (49) with the rest of the system residuals; this implies thai simulator and Hash residuals converge together without any special flash calculations. The other possibility, known as "exact flash" is to 10 determine equilibrium phases and compositions satisfying (49) to within a rigorous tolerance. In the exact flash approach, the residuals {49} are (numerically) zero due to flash iterations performed. Only derivatives of (49) with respect to simulator variables are required, to account for the contribution of equilibrium to the simulator Jacobian,

I S The exact flash policy is preferred, partly because it leads to a modular design, but principally because convergence of the highly nonlinear flash constraints frequently requires special treatment. In addition, the extra work inherent in performing exact Hash calculations is more than offsets by faster convergence of the simulator nonlinear iteration, resulting from more accurate equilibrium phases, 0

Efficient application of the stability testing and phase-split algorithms presented in sections 3,1 and 4 to the repeated solution of flash problems encountered in the course of simulator time-stepping requires that careful use be made of pre-existing conditions of pressure, temperature and composition on the simulation grid, 5

5.1 Two Phase Region

Generally speaking, saturated cells are overwhelmingly more likely to remain two- phase than to transition into single-phase in subsequent iterations. Furthermore, in the 30 vast majority of cases, the equilibrium phase compositions and -amounts already calculated provide an excellent starting point for new phase-split calculations at the updated cell conditions. Consequently, an attempt is always made to proceed directly to the Newu.ni procedure described in section. 4.3. using previous K-vaiυes as the initial guess and converging the governing equations (42 j to a strict residual tolerance {default isε_v,_ή, ~ 10^{' :}" ) in a very small number of iterations. Excessive iterations, or overshoot in the iteration variables, m the uπiikdy event that they should occur, results in swift termination of the iteration in favor o(Λb frrifio flash calculations.

5.2 Single Phase Region

^('ells which are undersaturattd are far more likely to remain single-phase than to evolve additional phases; however, this can only be established with complete certainty through stability testing, and the traditional approach requires this testing to commence from a crude starting point, such as the Wilson K- values, equation (29), using both a 'light' and a 'heavy^* trial phase, for each new condition encountered,

A worthwhile objective is to reduce the number of costly stability tests performed. This has been attempted in the past by limiting testing to certain candidate cells, such as those which border on clusters of existing two phase cells. However, such an algorithm is ultimately heuristic, and introduces a recursive component in the sense that (.nice instability has been revealed, new neighbors must be tested, This has undesirable implications in the parallel simulator.

The preferred approach which is related to that of Claus P, Rasmussen, KπstJan Krejbjerg, Michael S.,. Mieheisen and Kerϋti E. Bjurstrom, Increasing the Computational Speed of Flash Calculations with Applications for Compositional. Transient Simulations, Society of Petroleum Engineers, SPE 84181. February 2006 SPH Reservoir Evaluation & Engineering is adopted, which will be referred to here as Conditional Stability Tasting (CS T). Instead of considering loeai grid conditions, this approach is based on monitoring each single-phase cell's approximate location in the phase plane. The picture be low illustrates the cases that arise,

The singte-phase region \% subdivided into a "'shadow" region (C) next to the phase boundary, characterized as the set off and T for which the FPD equations Q¹S) have one non-trivia] solution with a corresponding positive value of the TPD. Beyond the shadow region is the remote region (D), in which only a trivial solution*^ ~ - ) to the TFD equations exists. The shadow zone thus acts as a buffer between states far into the single-phase region, and the two-phase reginn itself.

The central idea of CST is to attempt, to skip stability calculations in Λ>ne (D), provided the magnitude of change experienced in a given cell is sufficiently small. For cells in region (C) stability testing is "single-sided" and commences with Newton iteration from the previously calculated, nonuivial solution with positive TPD,

However, as the width of the shadow zone can be shewn to shrink to zero in the critical region, a measure of distance to the eritical point is needed in order \o safely bk^'φ calculations in zone (D). This measure, introduced by Michelsen. is given by the smallest Eigen value of the matrix

When conditions P\ T\ -] in a cell are found to lie in region (D), the matrix (SO) is calculated and its smallest Ei gen value Al = min(efg [ B)) determined. All values marked (*) are then stored for re-υse.

In subsequent calculations at new conditions /\ T, ς ₅ define state-variable changes with respect to base point as

ΔP - P - P^* , AT s 7 - 7\ Δz, - z, - z^* (5 ! ⁾

A conipittationai tolerance is established by sealing a user-supplied tolerance parameter (default value _<?< = 0.1 } by the tabulated, approximate distance from the critical point, έ: -^; £ .λl . The stability test can be skipped at new conditions, provided the following conditions are ali satisfied:

5.3 Ah lnitio Flash Calculations

In the procedures outlines in sections 5,1 and 5,2 above, die possibility for algorithmic failure exists. For example, an assumed two-phase state may in reality be single phase; or, the quality of ihe initial guess insufficient to allow convergence in Newton's method. Analogously, a single-phase state previously located in the shadow or remote region may have experienced a change in conditions which k too large to allow the inference that the fluid remains stable as a single phase. In such cases, arid far situations when no initial information is available, the Ab lnitio (lat, "from the beginning") flash calculation Is required. The present approach is based on Michelsen's work. The main difference relates to the use of the reduced variable technique of section 3.1 for the stability testing step.

I o begin the exposition, note that for two phases of assumed compositions x,>\ the difference in GFR with the fluid when viewed as single phase can be expressed as:

ΔG Ml - /O][X (log y, + log^ ) ± β∑y_. Oog^ *^■ ?og$^" } ^•■■ - ^r₁ (log ^, ÷ log#;^{" )}

Using raass-balance r, ~ (1 - β)x, + βy, this can be written

άG _^- O - β^ x_Λloi' X, + k>g>£^; - log s, " log^ κ /?]T v_; (k>g v -s- iog$^r - logs, ~ k>g$/^'')

or, recognizing the tangent -plane distances of the liquid and vapor compositions, ΔG - (1 - /?}TPDfx) + /?TPD(y)

If compositions x, v can be found, satisfying mass-balance, with ΔG < O _s the mixture ;:, is unstable at the current conditions .P and T . In addition, if 5 either TFD(x) < 0 or IPD(y) < 0the mixture is also unstable {in practice, a small threshold value is used instead of zero; default setting i$TPD < £_J(^ ~ -lG^{"! i} ). Tn such situations, stability testing is redundant.

The steps of the Ab InMo algorithm are as follows: K)

First, evaluate d. ( P, T, 2) , the CHbbs Free Energy of the feed and the Wilson K- values K " using equation (29). Next, perforin three cycles of successive substitution. if during this process ΔG < 0.0 is observed the feed is unstable ami the current compositions can he used in subsequent split calculations. If instead it is the case

15 that TPIXx) < Q or TPD(y) < Othe same conclusion applies. The selection of estimates will now be based on which phase indicated instability, e.g., IfTPD(X) -C O then estimate K-vaiues as log K₁ ^- log φ!' (x)~- log φ'^' . if after three iterations instability has not been revealed, no definite conclusion is possible and a full stability test is performed. 0 if instability is detected, either through stability testing or the SS approach outlined above, estimates are now available for which the objective function ΛG is negative. Three additional cycles, each consisting of three SS negations are applied, each followed by an attempt to accelerate the process. Only convergence tα a minimum of 5 the GFE is possible in this approach. If convergence tolerance is met following completion of the cycles, the algorithm terminates; if not, a 2ⁿ" order rigorous GFE minimization algorithm, enforcing strict descent, is applied for the final convergence. Consequently, the trivial solution is always avoided, and convergence to a minimum of the GFF; is guaranteed, 0

- 2& 5.4 Phase Labeling

Once it has been verified that the fluid of composition s, is stable as a single pha;>e at the prevailing pressure and temperature . a label of either 'oil⁵ or 'gas^7' must be 5 assigned to it, the principal reason being the need to apply the correct relative permeability table in calculating Slow properties of the phase,

ϊt is important to point out, however, that such a distinction gradually becomes meaningless as the critical point is approached, and that at present, no universally 10 accepted methodology for labeling exists - current simulators show a great deal of variability in this regard. From a physical standpoint, flow properties of a phase cannot reasonably depend on this label as the phases become indistinguishable.

While it seems clear that the most rigorous- approach to the labeling problem is the i -^■■ determination of the mixture true critical point this is a costly process, and an accurate determination may not be required, particularly if properties are extrapolated between oil and gas sear she true or approximate critical point.

For these reasons, the present simulator uses a simple correlation for pseudo -critical 20 temperature, svhich can be expected to be accurate enough to alimv correct labeling of the phase weil away from miscible conditions. This so-called Li-correlation represents a weighted average of the component, critical temperatures,

25 Here T is a correction factor that is typically unity, unless the model has het'O tuned to match initialization data, e.g. the location of a gas-oil contact.

Using (52), at tlie operating temperature 7^' a fluid satisfying T < T^₁₁ is labeled oil; otherwise, gas. Referring again So Figure 4 a method for reservoir simulation is illustrated, hi step

100, reservoiϊ model and reservoir data is imputed. The fluid phase properties and equations of state arc calculated in step ) 10. Such calculation;; can be those described in sections 1-4 herein. A Jacobian Matrix is generated in step 120 can be as taught in sections 1 -4.2 herein. The linear equations are then solved in step 130 pursuant to the teachings of sections 1-4,2.3, The solution is then updated in step S 40 pursuant to the teachings of section 1 -4,3. Then the solution is tested for stability or convergence pursuant to the teachings of sections 1 -4.2.3. Then the calculated soiuuun is output for the user in step 160.

Referring to Figure 5. a method 200 for reservoir simulation is illustrated, In step 210, a eel) is selected which has a vapor phase and a liquid phase within the cell. An estimated is made as to which of the vapor phase and the liquid phase is prcsem in a least abundant amount in step 220. In step 230.. the phase having the least abundant amount is assigned as the primary phase, and the other phase is assigned as the secondary phase. In step 240, the phase properties of the primary phase are computed utilizing the primary variables with the primary phase, and the phase properties υf the secondary phase arc also computed utilizing mass balance and die second variables associated with the secondary phase in step 250, stability is ensured in Use calculations by dividing by the primary phase rather then by a value near zero.

In the method, the phase properties of the priraaiy and secondary phases can include pressure, temperature, and pressure of the primary phase. In the method, the phase properties of ihe primary and secondary phases can include pressure, temperature, pressure, composition and amount of the primary phase.

in the method, the calculations of step 240 can farther cornpxise the following steps tor calculating the phase properties of the primary and secondary phases: (i) utilizing a reduced variable algorithm with the primary a;nd secondary variables associated with the primary and secondary phases, which produces a Rachford-Rice expression; (n) linearizing the Rach ford-Rice expression with K-values and reciprocal K~values. thereby creating linear expressions: (iii) generating a Jacυbian Matrix utilizing the primary and secondare variables; and (iv) solving the linear expressions and the Jacobiaπ Matrix to update the phase properties and test for stability. Such steps are taught in sections 1 -4,2.3 herein. In the method above, step t iii) can also include the steps oi^'; (I) calculating derivatives of the phase compositions; {2} calculating derivatives of the K-vaiues: and (?) calculating derivatives of fugacivy coefficients corresponding to each of the primary and secondary variables.

Referring to Figures 2 and 6, another method 300 is illustrated for determining the composition of fluid in a ceϊϊ during a computer enabled reservoir simulation, in step 310, direct reduced variable split calculations are performed using K- values when the cell had a fluid with a plurality of phases in a previous iiroeslep, In step 320, a single- sided reduced variable stability test is performed using vapor incipient moles when the ceil had o fluid with a single phase located in the shadow region, liquid side in the previous times tep. In step 33O₅ a single-sided reduced variable stability teat is performed using liquid incipient moles when the cell had a fluid with, a single phase located in the shadow region, vapor side m the previous timestep. In step 340. an ab MiIo Hash calculation is performed on the cell to determine fluid composition and proceed to step 360 when (he cell is in the remote region The ah htilks calculation can be as taught in section 5.3 herein. In step 350, there is a determination of whether there is a failure in steps 310, 320. or 330. In step 350, an ah srthio Hash calculation is also performed on the cell to determine fluid composition, and then the method proceeds to step 360 when it is determined there is a fail lire. Jn step 350, when the fluid is determined to be single phase, additional calculations are performed to determine location in phase plane, and steps 320-340 are repeated., In step 360, the calculated results are used when there is no failure.

Referring to Figure 7, within a system 390, a computer readable media 410 is illustrated that is utilized during a reservoir simulation for determining the composition of fluid in a cell. As will be readily appreciated by those skilled in the art, the computer readable media can also be a co.mpone«i of a system in which the computer readable media or software 410 interacts with an input device 400, such as a computer terminal, and a central processing unit (CPU) 4(50. Those skilled in the art will aim readily appreciate such a system can also be part of a computer network. The computer media 410 includes a data receiver 420 that receives input reservoir model and data from a source. The computer media 410 also includes a least abundant amount assignor 430 that estimates which of a vapor phase and a liquid phase of the fluid in the cell is present in a least abundaru amount responsive to the input data received by the data receiver. The least abundant amount assigner 430 also assigns the phase having the estimated 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 that computes phase properties of the primary phase utilizing the primary variables with the primary phase, The .Quid phase property calculator 440 also compxjtcs phase properties of the secondary phase utilizing mass balance and the second variables associated vvith the secondary phase. The fluid phase property calculator 440 thereby ensures stability by dividing by the primary phase rather then by a value near zero. The computer .media 410 also includes an output producer 450 that is adapted to produce and communicate the calculated phase properties of the fluid to a readable format, kn instance to a screen of a monitor or to a primer.

In the computer readable media 4K). the fluid phase property calculator 440 can also include s reduced variable algorithm .module 470 that produces a Raehford-Riee expression with the primary arid secondary variables associated with the primary and secondary phases. The Quid property calculator 440 can also include a linearizing module 480 that creates linear expressions from the Rach ford-Rice expression with K-vahses and redpiucal K-vaϊues. The fluid property calculator 440 can also include a jacobian Matrix generator 490 that generates a Jacobian Matrix utilizing the primary and secondary variables. 1 he fluid property calculator 440 can also include a solver and stability tester 500 that solves the linear expressions and the Jaeobian Matrix to update the phase properties and test for stability. The Jacobian Matrix generator 490 can also have a phase composition derivative submodule that calculates derivatives υl the phase compositions; a K-valm? suhmodnle that calculates derivatives of the K- ^■values; and a iugadty submodule that calculates derivatives of the ϊugacity coefficient corresponding to each of the primary and secondary variables.

Referring to Figures 2 and 6, another method 500 is illustrated for determining die composition of fluid in a cell during a computer enabled reservoir simulation. Step 510 is determining whether the cell had a single phase or a plurality of phases in a previous timestep. From step 510, other steps are performed depending upon where the determination in step 5 S O, Step 520 is performed if the eel! had a plurality of phases in the previous timestep. In step 520 direct reduced variable split calculations are performed using K-vaiues. Step 530 is performed if the cell had a single phase in the previous timestep and was located in the shadow region, liquid side in a phase plane. In step 530 a single-sided reduced variable stability test is performed using vapor incipient moles. Step 540 is performed if ihe ceil had a single phase in the previous time-step and was located in the shadow region, vapor side of the phase plane, in step 540, a single-sided reduced variable stability test is performed using liquid incipient moles. Step 560 is performed if the cell is in the remote region of the phase plane. In step 560 an ah wirio flash calculation on the cell is performed to determine fluid composition. In step 560, after performing the ab inirio calculation, the method then proceeds to step 580. Step 570 is performed after performing steps .520-540 in order to determine whether there is a failure in fast processing. If there is a failure, then an ah initio flash calculation is performed on the cell to determine fluid composition, if the fluid is found to he single phase, additional calculations are performed to determine location in phase plane to be used in subsequent iterations. Step SSO is then performed if there is no failure. In step 380 the calculated results are used.

While in the foregoing specification this invention has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, it will be apparent to those skilled in the art thai the invention is susceptible to alteration and that certain other details described herein can vary considerably without departing from she basic principles of the invention.

- .V} - For example, the prcseiu invention also includes a system and computer readable, media carrying instructions for performing a compositional reservoir simulation of a subterranean hydrocarbon-bearing reservoir. This system, including computer hardware and storage, wili carry out the method of reservoir simulation outlined 5 above. Similarly, the computer readable media carries instructions for performing a compositional reservoir simulation, of a subterranean hydrocarbon -bearing reservoir in accordance with the principles described above.

I O

Claims

^ΗΛX.JS.CLAlMlIliS:

L A methυU for reservoir simulation comprising;

(a) selecting a cell which has a. vapor phase and a liquid phase therewith! n; <b) estimating which of the vapor phase arid the liquid phase is present in a 5 least abundant amount;

(c) assigning the phase having the least abundant amoum as the primary phase and assigning the other phase as the secondary phase; (cl ) computing phase properties of the primary phase utilizing the primary variables with She primary phase, and i;onspuiiog phase properties of 10 the secondary phase utilizing mass balance and the second -variables associated with the secondary phase;

(e) ensuring stability by dividing by the primary phase father then by a value near zero.

2. The method of claim L wherein the phase properties of the primary phase ! 3 comprise pressure, temperature, and pressure of the primary phase,

3. The method of claim 1 , wherein the phase properties of the secondary phase comprise pressure, temperature, and pressure of the secondary phase,

4. The method of claim 1. wherein the phase properties of ihe primary phase comprise pressure, temperature, pressure, composition and amount oi the primary 0 phase.

5^, The method of claim 1 , wherein the phase properties of the secondary phase comprise pressure, temperature, pressure* composition and amounts of the secondary phase.

6. The method of claim 1 , wherein step fd) further comprises the following steps for calculating the phase properties of the primary and secondary phases:

(i) utilizing a reduced variable algorithm with the primary and secondary variables associated with the primary and secondary phases, which produces a Rachibrd-Rice expression:

(is) linearizing the Rachford-Rjce expression with K -values and reciprocal K-vaiues.. thereby creating linear expressions;

(iii) generating a Jacobian Matrix utilizing the primary and secondary variables: and

(iv) solving the linear expressions and the Jacobian Matrix to update the phase properties and test for stability.

7. The method of claim 6. wherein step (iU) further comprises:

{ 1 } eaioii! ating derivatives of the phase compositions; {2} calculating derivatives of the K- values; and (3) calculating derivatives of fugacUy coefficients corresponding to each of the primary and secondary variables. S. A method for determining the composition of fluid in a ceil during a computer enabled reservoir simulation, the method comprising the steps of;

(a) perform direct reduced variable split calculations using R values when the cell had a fluid with a plurality of phases in a previous timestep: (b) pes form a smgle-sided reduced variable stability test using vapor incipient moles when the eel! had a fluid with a single phase located in the shadow region, liquid side in the previous timestep; Cc) perform a single-sided reduced variable stability tet=t using liquid incipient .moles when the cell had a fluid with a single phase located in the shadow region, vapor side in the previous timestep;

(ά) perform an ab initio flash calculation on the cell to determine fluid composition and proceed to step (i) when the ceil is in the remote region;

(e) determine whether there is a failure in steps fa) (b) oτ fc^'i and φ perforπi an ub miiio flash calculation on the cell to determine fluid composition and pioceed to step (!) when it is determined there is a failure;

(i^'O perform additional calculations to determine location m phase plane and repeat steps (b) - (d;- when the fluid is determined to be single phase; and

(1) use ihc calculated results vvhen there is no failure.

9, A computer readable media that is utilised during a reservoir siinuiation ibr deterraiiiing the composition of fluid in a cell the readable media comprising; a data receiver that receives input reservoir model and data from a source;

,-> 7/ - a least abundant arøouut assignee that estimates which of a vapor phase and a liquid phase of the fluid in the ceil ss present in a least abundant araoum responsive 10 the input data received by the data receiver, and assigns the phase having the estimated least abundant amount as the primary phase and assigns the other phase as the secondary phase; a fluid phase property calculator that computes phase properties of the primary phase utilizing the primary variables with the primary phase, and computes phase, properties of the secondary phase utilizing mass balance and the second variables associated with the secondary phase, the fluid phase property calculator ensuring stability by dividing by the primary phase rather then by a value near aero; &.nά an output producer that is adapted to produce and communicate the calculated phase properties of the fluid to a readable format,

10 The computer readable media of claim % wherein the fluid phase property calculate r comprises ; a reduced variable algorithm module that produces a Rac^'hford-R.iee expression with the primary and secondary variables associated with the primary and secondary phases; a linearising module that creates linear expressions irorn the Rachford-Rice expression with K -values and reciprocal K -values; a Jaeobian Matrix generator that generates a Jacobian Matrix -utilizing ih« primary and secondary variables; and a solver and stability tester thaϊ solves the linear expressions and the Jacobian Matrix to update the phase properties and test for stability. 1 1 , The computer readable media of claim 10, wherein the Jacobiars Matrix generator further comprises: a phase composition derivative sishmodide that calculates derivatives of the phase compositions; 5 a fC-valae submocluie that calculates derivatives of the K-values; and a fαgacity submodule that calculates derivatives of the iugadty coefficients corresponding to each of the primary and secondary variables.

.12. A method for determining the composition of fluid in a cell di.iri.ng a computer enabled reservoir simulation, the method comprising the steps of-

} 0 (a) determining whether the ceil had a single phase or a plurality of phases in a previous time&tep; (ϊ) if the eel! had a plurality of phases m the previous timestep. then perform direct reduced variable split calculations using K- values;

1 5 (iij if the cell had a single phase in {he previous timestep, then take following action depending on that ceils previous location in a phase plane;

(1 ) i f the fluid is in the shadow region, liquid aide, ϋien perforin a single-sided reduced variable stability tesi 20 lifting vapor incipient moles;

{ 2) if the Quid is in the shadow region, vapor side, then perform a single-sided reduced variable stability test using liquid incipient moles. t iii) if She csli is in UK remote region, then perform an ab initio flash calculation on the cell to determine fluid composition and skip to step (c);

Ib) determine whether there is a faiiurc in fast processing by {i), (H), or (Hi);

Ci) If there is a failure, lhyti perform an ab hullo Hash calculation on the cell to determine fluid composition; (iij if the fiuid is found Jo be single phase, perform additional calculations to determine location in phase plane to be used in subseφitnl iterations;

(c) if there is no failure, then use the calculated results.