Generalized negative-flash method for multiphase multicomponent systems

Fluid Phase Equilibria 299 (2010) 272–284

Contents lists available at ScienceDirect

Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid

Generalized negative-flash method for multiphase multicomponent systems A. Iranshahr ∗ , D. Voskov, H.A. Tchelepi Department of Energy Resources Engineering, School of Earth Sciences, Stanford University, 367 Panama Street, Stanford, CA 94305, United States

a r t i c l e

i n f o

Article history: Received 6 February 2010 Received in revised form 13 September 2010 Accepted 14 September 2010 Available online 29 September 2010 Keywords: Negative flash Arbitrary number of phases Thermodynamic phase-state identification Tie-simplex Compositional space parameterization

a b s t r a c t We describe a generalization of the negative-flash method for computing the equilibrium compositions of systems that can form more than two fluid phases. We show the solution existence and uniqueness of the general Rachford–Rice problem for systems with any number of phases. A bisection based strategy is proposed either as a solution procedure to guarantee convergence, or as a preconditioning step of Newton based methods. A multi-stage negative-flash procedure is developed to identify the state of a multicomponent mixture where the maximum number of phases is arbitrary. The approach should be used in combination with a tie-simplex based modeling framework to provide a general method for phase behavior computations. Challenging numerical examples are presented to demonstrate the accuracy and robustness of the negative-flash approach. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Predictions of recovery by gas injection in oil reservoirs involve solving coupled systems of transport equations. Resolving the complex multiphase thermodynamic equilibrium behaviors associated with these subsurface displacement processes poses significant computational challenges [1–3]. Here, we propose a multiphase equilibrium method based on a generalized negative-flash strategy. Some compositional and many thermal-compositional displacement processes involve fluid mixtures that form three or more fluid phases at equilibrium. For such problems, the standard simulation approach based on two-phase equilibrium is inadequate. In compositional reservoir simulation, a stability test is performed to determine the phase-state of a mixture for a given overall composition at fixed pressure and temperature [4]. For multiphase states, a flash procedure is used to determine the phase amounts and compositions [5]. Flash calculations involve simultaneous solution of the thermodynamic constraints for all the components that may be present in the equilibrium phases. For gas injection processes, Method of Characteristics (MoC) solutions for two-phase multicomponent displacements have been derived for ideal systems (one-dimensional dispersion-free with constant pressure and temperature) [6]. In these MoC solutions, tie-lines are used to construct the solution path in compositional space. Motivated by this theory, a general tie-line based Compo-

∗ Corresponding author. E-mail addresses: [email protected] (A. Iranshahr), [email protected] (D. Voskov), [email protected] (H.A. Tchelepi). 0378-3812/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2010.09.022

sitional Space Parameterization (CSP) method has been developed in order to improve the computational efficiency and accuracy of isothermal compositional flow simulation [1,2]. The mathematical framework for the extension of the CSP methodology to an arbitrary number of phases has been presented in [3]. In two-phase displacements, compositions in the single-phase region connected to compositions in the two-phase region lie along key tie-lines that intersect the injection and initial compositions [1,2,6]. Tie-lines are used in the CSP approach to parameterize the compositional space in the vicinity of a mixture with a given overall composition [1,2]. Since the mixture can be in the single- or two-phase state, a two-phase negative-flash procedure is required. For displacements that involve more than two fluid phases, tie-simplexes, which represent mixtures at multiphase thermodynamic equilibrium, can be used to parameterize the compositional space [3]. In the generalized CSP approach, a robust multiphase negative-flash procedure is required, since the mixture can be outside the multiphase region enclosed by the parameterizing tie-simplex. The problem of how components partition among multiple phases at equilibrium was first formulated by Rachford and Rice [7]. Assuming constant partitioning coefficients (K-values), they derived an objective function to be solved for a phase fraction in a two-phase mixture. In their work, they only considered positive phase ratios. Later, Li and Nghiem [8] showed that the formulation can be applied to negative phase fractions, which allows phase split calculations for mixtures in the single-phase, as well as, the twophase region. They updated K-values in an iterative procedure to satisfy the equal fugacity constraints for each component. Whitson and Michelsen [9] made significant contributions to the negative-

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

flash problem. They bounded the range of the phase fractions for the Rachford–Rice problem, and they argued that application of two-phase negative-flash calculations for a single-phase mixture is as reliable as a traditional phase stability test. An advantage of the negative-flash procedure is that it calculates the tie-line that parameterizes a single-phase mixture. The negative-flash can also be used to determine the phase-state of a fluid mixture in a gridblock during flow simulation [9]. It also plays an important role in computing the Minimum Miscibility Pressure (MMP), a pressure value at which the given composition is on the extension of a critical tie-line [10]. Wang and Orr [11] derived a variant of the negative-flash for systems outside the compositional space (i.e., with negative compositions). Several authors have discussed thermodynamic stability testing for multiphase systems. Gautam and Seider [12] developed a phase splitting algorithm to test the stability of a mixture. Michelsen [4,5] introduced the tangent plane analysis for determining the thermodynamic stability of a fluid mixture. The problem of multiphase flash was tackled effectively by Michelsen [13]. In his work, Michelsen developed a procedure based on constrained minimization of a Gibbs energy function. The original method dealt with positive phase ratios, and it was later generalized to a negative-flash scheme by Leibovici and Nichita [14]. Here, we study the generalized Rachford–Rice problem for mixtures that can form an arbitrary number of phases. The algorithm combines Successive Substitution Iteration (SSI) with the Newton method [5]. During any SSI iteration, the K-values are assumed constant; therefore, a robust constant K-value negative-flash procedure (i.e., multiphase extension of the Rachford–Rice problem) is required to obtain a stable SSI algorithm. Our experience indicates that a gradient based methodology to solve the multiphase flash problem is not guaranteed to find the global minima [15]. The problem is, however, amenable to a bracketing based scheme, such as bisection. We show the existence and uniqueness of the solution for the multiphase negative-flash problem with constant K-values. Then, we provide a nested bisection solution procedure for three-phase systems. This procedure is different from the two-dimensional bisection method to solve the three-phase Rachford–Rice equations [16]. The method proposed here is a general negative-flash procedure extended recursively to systems with any number of phases. Accurate phase-state identification in systems that can form more than two fluid phases is challenging because any combination of phase-states is theoretically possible within the compositional space. Identification of the phase state using a multiphase negative-flash is sufficient only if the mixture is inside the tie-simplex. For three-phase systems, for example, performing a three-phase negative-flash is not enough to confirm the phase-state of single- or two-phase mixtures. Additional twophase negative-flash computations are required to identify the phase-state of such compositions. Our objective is to develop a negative-flash method that yields the correct phase-state as well as the compositions of the phases at equilibrium. One can show that tie-simplex compositions vary continuously with pressure, temperature, or along a continuous route in the compositional space. As a result, one can parameterize the tie-simplexes associated with a given feed at discrete pressure and temperature values [17]. Since the phase-state of the feed for which one tabulates the tie-simplexes is not known a priori, a robust negative-flash procedure is required. Given the continuity of the tie-simplex compositional space, interpolation as a function of pressure and temperature is used to determine the phase-state of a given composition. Theoretical analysis and extensive numerical simulations indicate that the number of tie-simplex tables constructed during a flow simulation is quite limited.

273

We demonstrate the robustness of the algorithm for a very wide range of physical parameters. Several examples of gas injection in oil reservoirs with complex multiphase behaviors are presented to demonstrate the effectiveness of the proposed method. 2. Multiphase negative-flash framework Here, we describe a general framework for solving the nonlinear set of equations in the multiphase negative-flash problem. Next, we review the two-phase problem, and formulate the negative-flash for systems with an arbitrary number of phases. 2.1. EoS based multiphase negative-flash We start with an overview of the multiphase negative-flash framework based on an Equation of State (EoS). Let Np represent the maximum number of the phases that can coexist for a given system at the specified temperature and pressure conditions. Here, we consider the partitioning of mixture made up of Nc components into Np phases at equilibrium. The set of the equations can be written as fˆi,0 − fˆi,j = 0,



i = 0, . . . , Nc − 1, j = 1, . . . , Np − 1,

(1)

Np −1

zi −

xi,j j = 0,

i = 0, . . . , Nc − 1,

(2)

j=0



Nc −1

(xi,0 − xi,j ) = 0,

j = 1, . . . , Np − 1,

(3)

i=0



Np −1

1−

j = 0,

(4)

j=0

where zi denotes the overall composition of component i. The fugacities in Eq. (1), fˆi,j , are computed using an EoS model (e.g., Peng–Robinson [18]). The above Np (Nc + 1) equations are solved for the fraction (j ) and composition (xi,j ) of each phase. The pure phases define the vertices of a tie-simplex in the compositional space. Noting that the actual number of phases at equilibrium may be smaller than Np , the phase fractions are not restricted to [0, 1]. Newton’s method can be used to solve Eqs. (1)–(4). Convergence of the Newton iterations depends on the initial guess. A good initial guess is obtained using a Successive Substitution Iterations (SSI) procedure [9]. During each iteration of SSI, component fugacities are held fixed, and Eqs. (2)–(4) are solved simultaneously. Here, we present a generalization of the two-phase negative-flash approach [9]. We note that a constrained optimization method can also be implemented [13,14]. Negative-flash iterations can be initiated using various techniques [3]. Here, we describe a general procedure to obtain an initial guess for the Nc -component Np -phase negative-flash problem. Consider an Mc -component system, where Np ≤ Mc ≤ Nc . The Np -phase negative-flash iterations can start from a tie-simplex that lies on a face of the Mc -component space. The problem of solving for a tie-simplex in the (Mc − 1)-component system can be tackled in a similar way, and the procedure is recursively extended to all lower dimensional compositional spaces with Mc ≥ Np . Performing an Np phase negative flash in an Np -component system is straightforward, since a single tie-simplex parameterizes the entire compositional space. Upon convergence, an Np -phase negative-flash provides the phase state of the mixture only if the fractions of all Np phases are positive (i.e., the mixture is in the Np -phase region). We show that all the Mp -phase regions with 1 ≤ Mp ≤ Np exist in the compositional space at given temperature and pressure conditions. Lower-dimensional tie-simplexes parameterize the compositional

274

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284 A 50 40

physical domain

30 20

a

10 f(ν1)

x y z

b

0

z

a

y

x

0

1

-10

b

-20 -30 -40 -50 -5

-4

-3

-2 ν

B

C

-1

1

a b

Fig. 1. (a) Two-phase region inside the compositional space of a three-component system with constant K-values and (b) f(1 ) in the physical domain of the problem. All the mixtures between a − x and y − b are thermodynamically stable (single-phase). For this system, K = {1.2, 1.8, 0.2}, and z = {0.45, 0.1, 0.45}.

space around a given tie-simplex. As a result, in order to identify the state of the compositions that are outside the Np -phase region, a multi-stage negative-flash framework can be developed. In general, an (Np − 1)-stage negative-flash is required for a system that can form a maximum of Np phases at equilibrium. At the ith stage, a set of (Np + 1 − i)-phase negative-flash calculations is performed, where the initial guess is obtained from a face of the calculated tiesimplexes at stage i − 1. The calculations are terminated once it is found that the mixture is inside a tie-simplex. The composition is single-phase if it is outside all the tie-simplexes calculated in the multi-stage procedure.

Noting that, ∂mi = gi = Ki − 1, ∂1

i = 0, . . . , Nc − 1,

(8)

the physical domain is a bounded set of phase fractions, 1 , if the K-values satisfy the following condition,

∃i, j such that gi gj < 0.

(9)

One can show that f(1 ) is a monotonically increasing function of 1 , and that there is only one unique solution. 2.3. Multiphase negative-flash formulation

2.2. Overview of the two-phase negative-flash In this section, we briefly review the two-phase negative-flash problem. Many of the ideas discussed here generalize to phase equilibrium calculations with more than two phases. Fig. 1(a) shows the phase diagram of a three-component mixture with constant K-values. Every composition on the bubble point locus (x) is at equilibrium with only one composition on the dew point locus (y). A tie-line connects all the compositions that share these two equilibrium mixtures. The same tie-line parameterizes all the single-phase compositions that lie on its extension (i.e., all compositions between a and b in Fig. 1(a)). By performing a negative-flash, a tie-line is computed for a given overall composition (z), regardless of its location in the single- or two-phase region. Allowing for negative equilibrium phase fractions with 0 + 1 = 1 [9], and expressing the overall compositions in terms of two equilibrium phase compositions, we write: zi = 0 xi,0 + 1 xi,1 = xi,0 + 1 (xi,1 − xi,0 ) = xi,0 mi (1 ), i = 0, . . . , Nc − 1,



f (1 ) =

i=0

 z (1 − K ) i

i=0

i

mi (1 )

= 0.

(6)

i = 0, . . . , Nc − 1,

i = 0, . . . , Nc − 1.

(10)

where we have arbitrarily chosen a phase as the ‘base’ phase, and we have denoted its index by 0. Using



Np −1

0 = 1 −

j ,

(11)

j=1

xi,j = Ki,j xi,0 ,

i = 0, . . . , Nc − 1; j = 1, . . . , Np − 1,

one obtains

⎛

zi = xi,0 ⎝1 +



Np −1

(12)

⎞ j (Ki,j − 1)⎠ = xi,0 mi (),

i = 0, . . . , Nc − 1. (13)

In order to determine the equilibrium partitioning of all components among all phases, one must solve the following coupled Np − 1 equations



Nc −1

fj () =

i=0

 zi (1 − Ki,j )

Nc −1

(xi,0 − xi,j ) =

i=0

mi ()

= 0,

j = 1, . . . , Np − 1, (14)

Fig. 1(b) shows f(1 ) in the space of physical phase fractions. From Eq. (5), the physical domain of the problem is defined by mi (1 ) > 0,

xi,j j ,

j=1

(5)

Nc −1

(xi,0 − xi,1 ) =



Np −1

zi = xi,0 0 +

j=1

where mi (1 ) = 1 + 1 (Ki − 1), and Ki = xi,1 /xi,0 . With constant Kvalues (the Rachford–Rice problem [7]), the objective function to be solved for the single unknown phase fraction, 1 , is Nc −1

We formulate the negative-flash procedure for multiphase systems with constant K-values. The overall mole fraction of a component, zi , is given by:

(7)

for Np − 1 unknowns, ( ≡ {j }, j = 1, . . ., Np − 1). Next, we define the domain of the full problem. For this purpose, let us assume that we have an admissible set of K-values (so that the domain of Eq. (14) is bounded). From Eq. (13), we conclude

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

275

A 4 3.5 3 boundary of physical domain (Ω)

2.5

ν2

2

D

1.5

3-phase region (Δ)

1

B

solution

0.5 0 -0.5 -2

C

-1

0

a

1

ν1

2

3

b

Fig. 2. (a) Three-phase region (wedge-shaped) inside the compositional space of a four-component system with constant K-values and (b) space of the physical phase fractions (˝), and the solution point associated with an arbitrary composition in (a). In this example, z = {0.412, 0.155, 0.369, 0.063}, and K = {0.938, 1.446, 0.543, 1.380;0.713, 0.953, 0.549, 3.533}.

that mi () must be strictly positive for each component (denoted by index i). Let

for the first Mp phase fractions assuming that the last Np − Mp − 1 phase fractions are known:

˝ = { ∈ RNp −1 |0 < mi (), i = 0, . . . , Nc − 1},

min j,˝ < j < max j,˝ ,

(15)

represent the set of ‘physical’ phase fractions. Given an admissible set of K-values, ˝ is a full-dimensional bounded polytope in RNp −1 [19]. Let  = { ∈ RNp −1 |0 ≤ j ≤ 1, j = 1, . . . , Np − 1},

(16)

define the bounded polytope representing the region where all Np phases coexist. We note that mi (ej ) = Ki,j > 0,

i = 0, . . . , Nc − 1, j = 1, . . . , Np − 1,

(17)

where ej is a unit-vector in RNp −1 . As a result,  ⊂ ˝. We next show that ˝ is a convex polytope in RNp −1 . Note that a line segment connecting any two arbitrary points in ˝ is itself within ˝. To see this, consider two points, 1 and 2 , inside ˝; and let 3 be a point on the line connecting 1 and 2 : 3 = ˛1 + (1 − ˛)2 ,

(18)



3,j (Ki,j − 1)

j



= 1+˛



1,j (Ki,j − 1) + (1 − ˛)

j

= ˛mi (1 ) + (1 − ˛)mi (2 ) > 0,

(20)

where j is a value in the range of the j-th phase fraction for the full-dimensional ˝. The equations to be solved simultaneously for ( ≡ {j }, j = 1, . . ., Mp ) are:

 zi (1 − Ki,j )

Nc −1

fj () =

i=0

= 0,

mi () + ni

j = 1, . . . , Mp ,

(21)

where ni is a fixed number given by, Np −1

ni =



j (Ki,j − 1),

i = 0, . . . , Nc − 1.

(22)

j=Mp +1

We note that mi () is a function of the Mp phase fractions only. Let ˝Mp = { ∈ RMp |mi () + ni > 0, i = 0, . . . , Nc − 1},

(23)

represent the set of physical phase fractions for this subproblem. If ˝ is bounded, ˝Mp will be bounded. Moreover, ˝Mp is convex. In order to show this, pick two arbitrary points, 1 , 2 ∈ ˝Mp , and a point on the line segment between them: 3 = ˛1 + (1 − ˛)2 , where 0 < ˛ < 1. Since

where 0 < ˛ < 1. Since mi (3 ) = 1 +

j = Mp + 1, . . . , Np − 1,

2,j (Ki,j − 1) + ˛ − ˛

mi (3 ) + ni = 1 +



j

i = 0, . . . , Nc − 1,

3,j (Ki,j − 1) + ni

j

(19)

3 ∈ ˝; and ˝ is convex. 3. Physical domain of the problem In this section, we show that the general negative-flash problem is recursive, and we describe the physical domain of the equations for a subproblem. Then, we specialize our developments for threephase systems.



= 1+˛



1,j (Ki,j − 1) + (1 − ˛)

j

2,j (Ki,j − 1)

j

+ ˛ − ˛ + ˛ni + (1 − ˛)ni = ˛ (mi (1 ) + ni ) + (1 − ˛) (mi (2 ) + ni ) > 0,

(24)

˝Mp

for i = 0, . . ., Nc − 1, is convex. Therefore, the subproblem of solving for Mp phase fractions has all the properties of a full negative-flash problem, and a recursive procedure can be applied to solve the multiphase negative-flash problem.

3.1. Solution domain for the recursive subproblem 3.2. Solution domain for three-phase systems We show the existence and uniqueness of the general negativeflash problem by considering subproblems of the full problem. Suppose for an Np -phase negative-flash problem, we want to solve

For three-phase systems, ˝ is a convex polygon, and we solve for two phase fractions (1 , 2 ). Fig. 2 describes the construction

276

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284 2

3

1.5 +

m (ν)=1+0.3ν -0.7ν

1

2 0.5

Δ

0

1

Ω

ν2

ν2

Ω

-0.5

+ m (ν)=1-0.1ν

ν =0

Δ

+

0

-1 m (ν)=1+0.2ν +0.5ν =0

-1.5

-1

-2

-2 -2.5 -5

-4

-3

-2 ν

-1

0

-2

1

-1

0

1 ν1

1

a

2

3

4

5

b

Fig. 3. ˝ for a 3-component system (a, unbounded), and a 25-component system (b, bounded) with arbitrary K-values. Note that in both examples,  ⊂ ˝.

of ˝ for a four-component system with an admissible set of Kvalues. The three-phase behavior inside the compositional space is shown using tie-triangles (see Fig 2(a)). We choose an arbitrary composition and find the unique tie-triangle that parameterizes this composition (i.e., its extension intersects the composition). We extend this tie-triangle plane throughout the compositional space to include all the compositions that it parameterizes. The bounds of this plane are shown by the dashed blue lines in Fig. 2(a). Since each composition on this plane is associated with a unique set of phase fractions, the set of independent phase fractions (1 and 2 ) is mapped into the phase-ratio space (˝). Note that ˝ is a convex bounded polygon (see Fig. 2(b)). In order to define the conditions that must be satisfied by the K-values, consider the following set of K-values for a threecomponent system,

 K=

0.9 2 1.2 1.5 1.3 0.3

 .

Fig. 3(a) shows the three lines where mi () = 0. Note that mi is positive on one side of each line (represented by the arrows). Although for this example  ⊂ ˝, ˝ is obviously not bounded. However, the setup suggests a procedure to verify the validity of a given set of K-values. We choose a point inside ˝,  = 0, and calculate the directional derivatives of mi () along Nc − 1 lines through 0 that are parallel to each mj =/ i () = 0. The equation of such a line is (Kj,1 − 1)1 + (Kj,2 − 1)2 = 0, and the phase fractions along this line change as follows



∂ ∂s



=



1



2

(Kj,1 − 1) + (Kj,2 − 1)

j

2

Kj,2 − 1 1 − Kj,1



,

(25)

with respect to a parameter s. Moreover,

dm () i

ds

j

= gi,j =



2  ∂m () ∂ i

∂k

k=1

=

k

∂s

2

(Kj,1 − 1) + (Kj,2 − 1)

,

(26)

2

determines how mi () varies when another mj () is fixed. Here, i and j refer to the component indices, i, j = 0, . . ., Nc − 1. For an unbounded ˝, there exists an mj () along which all other mi ()’s are either increasing or decreasing. Therefore, we deduce that ˝ is bounded if and only if

∃j, k such that gj,i gk,i < 0,

4. Existence and uniqueness of solution for the multiphase negative-flash problem Here, we prove that a solution to the set of equations exists within the physical domain, ˝, and that it is unique. We show this for the general multiphase case, and then we describe the developments for three-phase systems. 4.1. Existence and uniqueness of solution for the recursive subproblem We prove that by fixing some phase fractions (j with j = Mp + 1, . . ., Np − 1), the negative-flash problem for the rest of the phase ratios (j with j = 1, . . ., Mp , Eq. (21)) has one and only one solution. Let us first show the existence of such a solution. Since ∂fj ∂j

 zi (1 − Ki,j )2

Nc −1

=

i=0

(mi () + ni )2

i = 0, . . . , Nc − 1.

(27)

> 0,

j = 1, . . . , Mp ,

(28)

and given that fj () is indefinite at the boundaries of ˝Mp (Eq. (23)), the following equation has one and only one solution fj () = 0,

with fixed k =/ j ,

j = 1, . . . , Mp ,

(29)

where the subscript k = / j indicates that all the phase ratios, k , except j are held constant. As a result, fj () = 0 defines a surface which is extended over the range of j in ˝Mp . Since this is true for all of the equations (j = 1, . . ., Mp ), the surfaces defined by fj () = 0 must intersect for at least one point, and a solution to Eq. (21) exists. We next show the uniqueness of the solution. The Jacobian of Eq. (21) is

j

(Ki,1 − 1)(Kj,2 − 1) − (Ki,2 − 1)(Kj,1 − 1)



In this procedure, the number of constraints which must be honored is the number of components. Fig. 3(b) shows an example of bounded ˝ for a 25-component system with arbitrary K-values.

Ji,j =

∂fi = ∂j

 zk (1 − Kk,i )(1 − Kk,j )

Nc −1

k=0

(mk () + nk )2

,

i, j = 1, . . . , Mp .

(30)

We note that J is a symmetric matrix. First, we show that J is positive definite. For this purpose, we need to prove: 1. uT Ju ≥ 0, ∀ u ∈ ˝Mp , 2. uT Ju = 0 implies u = 0. We have,

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284 Mp Mp  

uT Ju =

277

2

uj uk Jj,k

j=1 k=1

zi

Mp 

=

i=0

(mi () + ni )2

j=1



Nc −1

=

zi p2i (mi () + ni )2

i=0

(1 − Ki,j )uj

Mp 

(1 − Ki,k )uk

solution

1

k=1

,

f2(ν)=0

0.5

(31)

0

where, pi =

f1(ν)=0

1.5

ν2



Nc −1

Mp 

(1 − Ki,j )uj .

-0.5

(32)

j=1

-1.5

-1

uT Ju ≥ 0.

Moreover, since both zi and mi () + ni are Therefore, strictly positive, uT Ju = 0 implies that p =0. From Eq. (32), this is the case if and only if u =0. As a result, J is a symmetric positive definite matrix. In order to show the uniqueness of the solution, we study how the value of fMp () changes along the curve defined by fj () = 0,

j = 1, . . . , Mp − 1,

(33)

in ˝Mp . In order to prove the uniqueness of the solution, it is necessary and sufficient to show that the directional derivative of fMp () with respect to Mp is strictly positive along this trajectory. Differentiating Eq. (33) yields:



Mp −1

∂fj ∂Mp

=−

Jj,k

k=1

j = 1, . . . , Mp − 1.

dMp

 ∂fMp ∂ k

Mp −1

=

+

∂k ∂Mp

k=1

JMp −1

J = JMp = where Mp −1

J





=

∂fMp ∂Mp

,

(35)



,

(36)

∂Mp

,

∂k



,

∂Mp

0.5

positive. Since the numbering of the first Mp phases is arbitrary, dfj dj

> 0,

along fk =/ j = 0,

k, j = 1, . . . , Mp ,

For the three-phase negative-flash problem, mi () = 0 is represented by a line in the space of phase fractions. Each line is associated with one component. The corners of ˝ are determined from the intersection of such lines. Consider the intersection of the pth and qth lines. This intersection is considered to be a corner of ˝, if it satisfies:

j, k = 1, . . . , Mp − 1,

(41)

for all the component indices (i = 0, . . ., Nc − 1) except p and q. Once all the corners of ˝ in the phase-ratio space have been computed, the minimum and maximum values of 1 and 2 can be determined (˝ is a convex polygon). For a given 1 such that 1,min < 1 < 1,max , one can calculate the range of physical 2 . Let j and k represent

From the Sylvester’s criterion [20], both JMp and JMp −1 have positive determinants. Noting that JMp −1 is invertible, Eq. (34) can be written as −JMp −1 w = v,

(37)

 Mp −1 −1

w = −(J

)

v=



∂j ∂Mp

,

j = 1, . . . , Mp − 1,

(38)

defines how phase ratios change along the trajectory. Consequently, from Eq. (35), dfMp dMp

 ∂f ∂ k k

Mp −1

=

k=1

∂Mp ∂Mp

(40)

and the solution to the subproblem is unique.

j = 1, . . . , Mp − 1.

where

1

Fig. 4. ˝ for a 15-component mixture with random K-values. Zero loci of both functions are shown by dashed curves. Their intersection is the solution, which lies outside .

mi () ≥ 0,



∂fj

∂fj

v=

v ∂fMp

vT

0

(34)

along the trajectory defined by Eq. (34). Let



ν1

4.2. Solution existence for three-phase systems

∂k , ∂Mp

We need to determine the value of dfMp

-0.5

+

= −vT (JMp −1 )−1 v +

∂fMp ∂Mp ∂fMp

∂Mp

= vT w +

=

∂fMp ∂Mp

det J Mp . det J Mp −1

(39)

We conclude from the Sylvester’s criterion that dfMp /dMp is strictly

Fig. 5. f1 () and f2 () at zero loci of the other function.

278

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

subsets of component indices for which Kj,2 > 1, and Kk,2 < 1, respectively. The physical bounds of 2 are obtained from,



max



1 (1 − Kj,1 ) − 1

≤ 2 ≤ min

Kj,2 − 1

j



k

1 (1 − Kk,1 ) − 1 Kk,2 − 1



.

Therefore,

N −1 2 c  zi (1 − Ki,1 )(1 − Ki,2 )

(42)

Similarly, the physical domain for 1 can be calculated for a given 2 ∈ (2,min , 2,max ). Fig. 4 shows the physical boundaries of the two phase fractions (˝, shown by black lines), as well as  for a 15component mixture with an arbitrary set of K-values that satisfy Eq. (27). From Eq. (14), we have: ∂fj () ∂j

 zi (1 − Ki,j )2

Nc −1

=

m2i

i=0

∂f2 () ∂f1 () = = ∂2 ∂1

,

j = 1, 2,

(43)

 zi (1 − Ki,1 )(1 − Ki,2 )

Nc −1

m2i

i=0

.

Here, we show that f1 = 0 and f2 = 0 intersect at one and only one point, i.e., Eq. (14) have a unique solution. For this purpose, we need to show that the directional derivative of either function along the zero locus of another function is strictly positive (or negative). We have df2 (1 , 2 ) ∂f2 (1 , 2 ) ∂f2 (1 , 2 ) ∂1 = + , d2 ∂2 ∂1 ∂2

(45)

where ∂ 1 /∂2 is determined by the locus of f1 (1 , 2 ) = 0. That is df1 (1 , 2 ) ∂f1 (1 , 2 ) ∂f1 (1 , 2 ) ∂1 = + = 0, d2 ∂2 ∂1 ∂2

(46)

(48)

where the derivatives on the right hand side are given by Eqs. (43) and (44). Note that Eq. (48) can be obtained from Eq. (39) with Np = 3 and Mp = 2. Consider the following two vectors in Nc -dimensional space: √ zi (1 − Ki,1 ) vi = , (49) mi √ zi (1 − Ki,2 ) wi = . (50) mi From the Cauchy–Schwarz inequality, ≤

i=0

i=0

 wi2

(52)

fˆi,0 = fˆi,j ,

i = 0, . . . , Nc − 1, j = 1, . . . , Np − 1.

(53)

where the fugacity of component i in phase j, fˆi,j , is a function of pressure, temperature, and the molar amounts of all the components present in phase j. In a negative-flash problem, we solve Eq. (53) for a mixture with fixed overall compositions, temperature, and pressure. Taking the molar amounts of the phases with j = 1, . . ., Np − 1 as independent variables, we differentiate Eq. (53) for the negative-flash problem to obtain:

 ∂ ln fˆi,0

Nc −1

∂nl,0

l=0

 ∂ ln fˆi,j

Nc −1

dnl,0 = p,T,nk = / l,0

∂nl,j

l=0

dnl,j ,

(54)

p,T,nk = / l,j

for i = 0, . . ., Nc − 1, j = 1, . . ., Np − 1. The subscript nk =/ l,j indicates that the molar amounts of all the components except l are held fixed in the j-th phase. Substituting



Np −1

dnl,m

(55)

ıj,m dnl,m

(56)

m=1

dnl,j =

 m=1

in Eq. (54) yields

2

v2i

∂f1 (1 , 2 ) ∂f2 (1 , 2 ) . ∂1 ∂2

Np −1

(∂f1 (1 , 2 )/∂2 ) , ∂f1 (1 , 2 )/∂1

N −1  N −1 c c  

≤

In this section, we develop a convergence criterion for the Successive Substitution Iterations of the multiphase negative-flash problem. This analysis is a generalization of the developments for the two-phase negative-flash problem [5,9]. Let us consider the equality of fugacities for each component between every pair of fluid phases:

and

df2 (1 , 2 ) (∂f2 (1 , 2 )/∂1 )(∂f1 (1 , 2 )/∂2 ) = ∂f2 (1 , 2 )/∂2 − d2 ∂f1 (1 , 2 )/∂1 = ∂f2 (1 , 2 )/∂2 −

2

,

m2i

i=0

Moreover, since vi and wi are linearly independent, the equality in Eq. (52) does not hold. Therefore, df2 (1 , 2 )/d2 (and similarly df1 (1 , 2 )/d1 ) is strictly positive (see Fig. 5), and the solution to Eq. (14) is unique.

(47)

Substituting Eq. (47) in Eq. (45) yields

vi wi

∂f1 (1 , 2 ) ∂2

dnl,0 = −

∂1 −(∂f1 (1 , 2 )/∂2 ) = . ∂2 ∂f1 (1 , 2 )/∂1

i=0



m2i

i=0

5. Convergence analysis

4.3. Solution uniqueness for three-phase systems

2

≤

N −1  N −1  c c   zi (1 − Ki,1 )2 zi (1 − Ki,2 )2

(44)

From Eq. (43), ∂ f1 /∂1 > 0 and ∂ f2 /∂2 > 0. Moreover, both f1 and f2 are indefinite at the boundaries of the compositional space. As a result, for a given 2 , there is one and only one 1 such that f1 (1 , 2 ) = 0. Similarly, for a given 1 , there is one and only one 2 , such that f2 (1 , 2 ) = 0. Therefore, the loci of f1 () = 0 and f2 () = 0 are represented by two curves extending from 2,min to 2,max and from 1,min to 1,max , respectively (see Fig. 4). These curves intersect for at least one point, which is a solution to Eq. (14). We note that the solution need not be inside the three-phase region ().

N −1 c 

m2i

i=0

   ∂ ln fˆi,0

(51)

∂nl,0

l=0 m=1

+ ıj,m

∂ ln fˆi,j ∂nl,j

p,T,nk = / l,0





dnl,m = 0, p,T,nk = / l,j

(57) where i = 0, . . ., Nc − 1, j = 1, . . ., Np − 1, and ıj,m is the Kronecker delta. Differentiating



Nc −1

ln fˆi,j = ln ni,j − ln

ni,j + ln ϕ ˆ i,j + ln p

(58)

i=0

yields



∂ ln fˆi,j ∂nl,j

.



Nc −1Np −1

p,T,nk = / l,j

ıi,l = − ni,j

1



Nc −1

i=0

ni,j

+

∂ ln ϕ ˆ i,j ∂nl,j

(59) p,T,nk = / l,j

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

279

1

C

1

three-phase

0.9

single- or two-phase

0.8 0.7

λ

0.6 0.5 0.4 0.3

p=211.05 bar

0.2 p=75 bar

0.1 0 NC

H2O

10

80

100

120

a

140 p

160

180

200

b

C

1

1

0.9

three-phase

single- or two-phase

0.8 0.7 0.6 λ

T=500 K

0.5 0.4 0.3 0.2 0.1

NC

H2O

10

0 500

510

520

530

540

T

T=547.9 K

c

d

Fig. 6. Continuous variation of the parameterizing tie-triangles for {C1 (0.2), NC10 (0.4), H2 O(0.4)} at (a) T = 500 K, and (c) p = 75 bar until degeneration. The largest eigenvalues of M for (b) constant temperature, and (d) constant pressure approach unity with increase in p and T, respectively.

As a result, the Jacobian of Eq. (57) is obtained from E = Ji,j;l,m

ıi,l − ni,0

1



Nc −1



+

∂ ln ϕ ˆ i,0

It can be shown that the error vector for Eq. (57) at iteration k + 1, ek+1 , is related to ek by



∂nl,0

ek+1 = Mek = (JA )

p,T,nk = / l,0

ni,0

i=0

⎡

⎤

⎢ ⎢ı 1 ⎢ i,l + ıj,m ⎢ − n ⎢ i,j N c −1 ⎣

1 + 2



∂ ln ϕ ˆ i,j



∂nl,j

ni,j

⎥ ⎥ ⎥ ⎥, ⎥ p,T,nk = / l,j ⎦

(60)

i=0

for i, l = 0, . . ., Nc − 1, and j, m = 1, . . ., Np − 1. Noting that during each Successive Substitution Iteration, the fugacity coefficients are composition-independent, an approximation of the exact Jacobian in Eq. (60) is obtained as

⎡

A Ji,j;l,m

ıi,l = − ni,0

1



Nc −1

i=0

ni,0

⎤

⎢ ⎢ı 1 ⎢ i,l + ıj,m ⎢ − ⎢ ni,j N c −1 ⎣

ni,j

i=0

⎥ ⎥ ⎥ ⎥. ⎥ ⎦

−1

(JA − JE )ek ;

(62)

as a result, a necessary condition for the convergence of the iterations is that all the eigenvalues of M must be smaller than one in magnitude [15]. We demonstrate this using a numerical example (see Fig. 6). We consider the three-phase negative-flash of a three-component mixture at fixed temperature and pressure. All the eigenvalues of M are calculated using the negative-flash solutions. In the first example, we plot the eigenvalues as a function of pressure. Note that two sides of the tie-triangle converge at p = 211.05 bar, and consequently, the two largest eigenvalues of M approach unity. However, for changes in temperature, the parameterizing tie-triangle degenerates (i.e., its area becomes zero) with three distinct phases. This is confirmed by the fact that all the eigenvalues of M remain smaller than one. 6. Numerical implementation and results

(61) In this section, we describe the implementation of the negativeflash procedure for three-phase problems. Then, we present a set of numerical results for three- and four-phase systems.

280

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

∞

∞

2

2

1.5

1.5

1

0.27852

1 2

-0.089836

ν

ν

2

-0.089836 0.5

0.5

0

0

-0.5

-0.5

-∞ -1.5

-1

-0.5

ν

0

0.5

1

-1.5

-1

-0.5

1

it er . 1

ν

0

0.5

1

0

0.5

1

1

i te r . 2

2

2

0.27852

1.5

1.5

0.027717

1

2

ν

ν

2

-0.089836 0.5

0.5

0

0

-0.5

-0.5 -1.5

-1

-0.5

0.027717 -0.032952 -0.089836

1

ν

0

0.5

1

1

it er . 3

-1.5

-1

-0.5

ν

1

i te r . 4

Fig. 7. The first four iterations of the methodology for the system shown in Fig. 4. The numbers represent the values of f2 (1 , 2 ) at 2,min , 2,mid , and 2,max where f1 (1 , 2 ) = 0. Note that initially the value of f2 is indefinite at 2,min and 2,max (see Fig. 5). These lower and upper bounds become closer as the iterations proceed.

6.1. Bisection based algorithm for three-phase systems Based on the properties of the problem discussed previously, a bisection based method is guaranteed to converge to the solution. In order to improve the speed of convergence, the bisection strategy can be combined with a gradient based scheme. Here, we use a nested bisection scheme to solve for the two unknowns 1 , and 2 in a three-phase negative flash problem. This algorithm is summarized as follows: Step 1: Given the set of K-values, find the boundaries of the physical plane, which yield 2,min and 2,max . Step 2: Let 2,mid = (2,min + 2,max )/2. Step 3: Solve f1 (1 , 2,min ) = 0, and f1 (1 , 2,mid ) = 0 for 1,min , and 1,mid , respectively using bisection. Step 4: If f2 (1,min , 2,min ) × f2 (1,mid , 2,mid ) > 0, 2,min = 2,mid ; else 2,max = 2,mid . Step 5: Stop if 2,max − 2,min < ε; otherwise, go to Step 2. Fig. 7 illustrates the first iterations of the procedure for the 15component system introduced previously (Figs. 4 and 5). The SSI method for solving two-phase negative-flash problems can be generalized to three-phase systems by introducing this algorithm as a constant K-value flash solver at each iteration. An initial guess for the K-values is required for such an SSI procedure. For this purpose, we flash the overall composition which yields the ‘largest’ tie-triangle at the given temperature and pressure conditions. The resulting tie-triangle is then used as an initial guess to flash the tar-

get overall composition. This method of computing an initial guess is applicable for systems with an arbitrary number of phases [3]. 6.2. Multi-stage three-phase negative-flash The three-phase negative-flash introduced in the previous section is not sufficient for determining the phase-state of a given composition. This is because the actual two-phase behavior outside the three-phase region may be different from the three-phase negative-flash results (the extension of the edges that define the tie-triangle, see Fig. 8). Based on a three-phase negative-flash only, a mixture can be identified as stable in the single-phase region; however, its Gibbs free energy might decrease by splitting into two mixtures. To avoid such problems, a second stage using two-phase negative-flash can be employed. Consequently, in order to determine the number and compositions of the existing phases for a given overall composition at specified temperature and pressure conditions, we use a multistage negative-flash algorithm. In the first stage, a three-phase negative-flash is performed. If the composition is inside the tietriangle (three-phase), its state is confirmed. Otherwise, based on the results of the three-phase negative-flash, additional two-phase negative-flash computations are necessary. Here, we assume that the tie-lines are in the plane of the tie-triangle, even though real multicomponent systems show slight deviations from this behavior. The compositions of the tie-triangle edges provide good initial guesses for the two-phase negative-flashes. This is because tieline compositions change continuously around a tie-triangle edge.

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

C2

C1

H2S

Fig. 8. Complete parameterization of {C1 , C2 , H2 S} at T = 180 K, p = 15 bar. The accurate two-phase behavior is shown by the solid curves outside the three-phase region.

281

A tie-line can eventually degenerate to a point (see Fig. 8). Tiesimplexes also degenerate with increase in pressure or temperature (see Fig. 9). We present two numerical examples to illustrate the parameterization of three-phase compositional spaces (see Fig. 10). For each system, we choose two compositions for the initial oil and injection gas. We identify the state and parameterize the tie-triangle/tie-line for the compositions along the dilution line (i.e., the line connecting the initial and injection compositions). One tie-triangle and two tie-lines are used to parameterize a single-phase composition. A two-phase composition is parameterized by one tie-triangle and one tie-line. Only a single tie-triangle is required to determine the state of a three-phase composition. In Fig. 10, phase boundaries show where the phase-state changes along the dilution line between the injection and initial compositions. The three-phase negative-flash procedure can also be used to generate phase diagrams (see Fig. 11), and this is demonstrated using two examples: (1) a sour gas example, and (2) a Maljamar reservoir oil [16]. The accuracy of our numerical results is verified by a commercial simulator [21] (For compositions and thermodynamic properties of these mixtures, refer to Appendix A). We consider dilution of both mixtures with pure CO2 . That is, for every mixture which results from the addition of CO2 to the initial

C1

C

1

C

H2O

10

C10

H2O

a

b

Fig. 9. Complete parameterization of {H2 O, C1 , C10 } at T = 520 K, (a) p = 100 bar and (b) tie-triangle degeneration pressure, p = 170.53 bar

CO2

H O 2

C1

C6 C

1

C10

C

CO2

2

a

b

Fig. 10. (a) Injection of CO2 into {C1 (0.1), C10 (0.4), H2 O(0.4), CO2 (0.1)} at T = 500 K and p = 140 bar and (b) injection of {C1 (0.5), CO2 (0.5)} into {C6 (0.5), CO2 (0.2), C2 (0.3)} at T = 180 K and p = 4 bar. Linear path includes compositions in gas-phase ( ), oil-phase ( ), two-phase (•), and three-phase ( ) states. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

282

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

a 40 b

100

35 95

30

90

p (bar)

p (bar)

25 20 15

85

80

10 75

5 70

0

0.3

0.4

0.5

0.6

0.7

0.8

0.65

0.7

0.75

z CO

0.8

0.85

0.95

0.9

1

z CO

2

2

Fig. 11. Phase diagram for the three-phase region (solid), and the locus of tie-triangle degeneration (dashed): (a) Sour gas example, T = 178.8 K and (b) Maljamar reservoir oil, T = 305.4 K.

feed, we increase the pressure and record the three-phase behavior. The pressure is increased until the parameterizing tie-triangle degenerates to a tie-line. Fig. 11 shows the locus of the threephase boundary, as well as, the tie-triangle degeneration locus. We note that no phase-stability test is required, since the threephase negative-flash is capable of determining the bounds of the three-phase region accurately.

H2O

6.3. Four-phase systems Given that the negative-flash subproblem has a unique solution, and that all the subproblems of a multiphase negative-flash problem have similar properties, we employ a recursive methodology. For a four-phase system, for example, we solve for 3 at the highest level (Mp = 3). During these iterations, once the values of 3 are obtained, the first two equations are solved (for 1 , 2 ) at the next (lower) level of the recursion (Mp = 2, see Eq. (21)). Similar to the three-phase negative-flash problem, one can derive the conditions for admissible K-values of four-phase systems. Here, ˝ is a three-dimensional convex polytope. In order to show that ˝ is bounded, one needs to prove that the Nc planes through  = 0 which are parallel to the mi () = 0 planes, are bounded. This is very similar to verifying the validity of a set of K-values for the three-phase negative-flash problem. Next, we present the compositional space parameterization approach for four-phase systems. Once the tie-tetrahedron, which encloses the four-phase region in the compositional space, is computed, we parameterize the lower-dimensional phase regions using tie-triangles and tie-lines. Fig. 12 shows an example of the fourphase behavior. We note that since in these conditions, water is

N2 C2

CO

2

Fig. 12. Four-phase compositional space parameterization at T = 140 K and p = 150 bar. The figure shows one tie-tetrahedron (in red), and three sets of side tie-triangles (in blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

almost immiscible with the other three phases, one side of the tietetrahedron almost lies on the C2 –N2 –CO2 face. Fig. 13 shows an example for which three-phase space parameterization alone is not sufficient to describe the compositional phase behavior [3].

CO2

CO2

N2

N2

C2 H S 2

C2

H2S

Fig. 13. Four-phase compositional space parameterization at T = 140 K and p = 200 bar. Each figure shows the tie-tetrahedron and two sets of side tie-triangles.

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284 Table A.1 Thermodynamic properties of the sour gas example.

20 18 16 14 12

p (bar)

283

Component z

Tc (K) Pc (bar) ω

Vc (L/mol) ıCO2

CO2 N2 C1 C2 C3 H2 S

304.2 126.2 190.6 305.4 369.8 373.2

0.094 0.0895 0.099 0.148 0.203 0.0985

0.2392 0.2334 0.3593 0.1011 0.067

73.8 33.5 45.4 48.2 41.9 89.4

0.225 0.04 0.008 0.098 0.152 0.081

ıN2

0.00 −0.02 −0.02 0.00 0.103 0.031 0.13 0.042 0.135 0.091 0.096 0.176

ıH2 S 0.096 0.176 0.080 0.070 0.070 0.000

10 Table A.2 Thermodynamic properties of the Maljamar reservoir oil.

8 6

Component

4

CO2 C1 C2 C3 NC4 C5–7 C8–10 C11–14 C15–20 C21–28 C29+

2 0

0.05

0.1

0.15

zCO

0.2

2

Fig. 14. Phase diagram for the four-phase region: sour gas example with T = 123.15 K.

We finally show how our negative-flash methodology can be used to construct four-phase diagrams (see Fig. 14). By decreasing the temperature of the sour gas example discussed before [16], four-phase behavior is revealed for small injection fractions. From this example, we note that by increasing the pressure, it is possible for a single composition to experience the three- and four-phase regions more than once. 7. Conclusions We formulated and solved the Rachford–Rice problem for an arbitrary number of fluid phases. We derived the conditions that a set of K-values must satisfy, and we proved the existence and uniqueness of the solution to the general problem. Based on these ideas, we developed a recursive bisection methodology to guarantee convergence of our algorithm. Due to its importance in reservoir simulation applications, we focused on the three-phase case. Our results indicate that the Np -phase negative-flash can be used to identify accurately the phase-state of mixtures that partition into Np equilibrium phases. For mixtures that partition into a lower number of phases at given pressure and temperature conditions, we proposed a multi-stage negative-flash procedure. Several fourphase examples are also presented to illustrate the methodology. List of symbols

z

Tc (K)

Pc (bar)

ω

Vc (L/mol)

ıCO2

0.2939 0.1019 0.0835 0.0331 0.1204 0.1581 0.0823 0.0528 0.0276 0.0464

304.2 190.6 305.4 369.8 425.2 516.7 590.0 668.6 745.8 812.7 914.9

73.8 45.4 48.2 41.9 37.5 28.8 23.7 18.6 14.8 12.0 8.5

0.225 0.008 0.098 0.152 0.193 0.265 0.364 0.499 0.661 0.877 1.279

0.094 0.099 0.148 0.203 0.255 0.3998 0.5373 0.7414 0.9815 1.2164 1.6267

0.000 0.103 0.130 0.135 0.130 0.103 0.103 0.103 0.103 0.103 0.103

Table A.3 Thermodynamic properties for C6 , NC10 , H2 O. Component

Tc (K)

Pc (bar)

ω

Vc (L/mol)

ıCO2

C6 NC10 H2 O

507.4 617.6 647.3

29.69 21.08 220.48

0.296 0.49 0.344

0.37 0.603 0.056

0.125 0.11 0.2

Table A.4 Binary interaction coefficients with H2 O. Component

ıH2 O

C1 C2 N2 NC10

0.4907 0.4911 0.275 0.45

ε J

tolerance for convergence of iterations Jacobian matrix

Acknowledgements We thank the Stanford University Petroleum Research Institute for Reservoir Simulation (SUPRI-B) for financial support. Appendix A. Thermodynamic properties of mixtures

Nc Np ˝ ˝Mp  j ıi,j fˆi,j ϕ ˆ i,j ni,j xi,j zi Ki,j 

number of components number of phases physical domain of the full negative-flash problem physical domain of a recursive subproblem Np -phase region in ˝ molar ratio of phase j the Kronecker delta, ıi,i = 1; ıi,j = 0, if i = / j fugacity of component i in phase j fugacity coefficient of component i in phase j molar amount of component i in phase j molar fraction of component i in phase j overall molar fraction of component i equilibrium factor for component i in phase j eigenvalue of M in Eq. (62)

The following four tables provide all the information required in the calculation of our numerical examples (Figs. 6–14). We have used the Peng–Robinson’s equation of state throughout our examples [18]. Tables A.1 and A.2 represent the molar compositions and properties of the sour gas and the Maljamar reservoir oil mixtures, respectively [16]. The binary interaction coefficients between hydrocarbon components [21] are calculated from:

 ıi,j = 1 −

1/6

1/6

2Vc,i Vc,j 1/3

1/3

Vc,i + Vc,j

1.2 .

For other thermodynamic properties, see Tables A.3 and A.4.

(A.1)

284

A. Iranshahr et al. / Fluid Phase Equilibria 299 (2010) 272–284

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

D. Voskov, H.A. Tchelepi, SPE J. 14 (3) (2009) 431–440. D. Voskov, H.A. Tchelepi, SPE J. 14 (3) (2009) 441–449. D. Voskov, H.A. Tchelepi, Fluid Phase Equilib. 283 (2009) 1–11. M.L. Michelsen, Fluid Phase Equilib. 9 (1982) 1–19. M.L. Michelsen, Fluid Phase Equilib. 9 (1982) 21–40. F.M. Orr Jr., Theory of Gas Injection Processes, Tie-Line Publications, Copenhagen, Denmark, 2007. H.H. Rachford, J.D. Rice, Pet. Trans. AIME 195 (1952) 237–238. Y.K. Li, L.X. Nghiem, SPE (1982) 11198. C.H. Whitson, M.L. Michelsen, Fluid Phase Equilib. 53 (1) (1989) 51–71. Y. Wang, Analytical calculation of minimum miscibility pressure, Ph.D. Dissertation, Stanford University, 1998. Y. Wang, F.M. Orr Jr., Fluid Phase Equilib. 139 (1) (1997) 101–124.

[12] [13] [14] [15] [16] [17]

[18] [19] [20] [21]

R. Gautam, W.D. Seider, AIChE J. 25 (1979) 999–1006. M.L. Michelsen, Comput. Chem. Eng. 18 (1994) 545–550. C.F. Leibovici, D.V. Nichita, Fluid Phase Equilib. 267 (2) (2008) 127–132. M.L. Michelsen, J.M. Mollerup, Thermodynamic Models: Fundamentals & Computational Aspects, Tie-Line Publications, Copenhagen, Denmark, 2004. K.B. Haugen, L. Sun, A. Firoozabadi, SPE (2007) 106045. A. Iranshahr, D. Voskov, H.A. Tchelepi, Theory and practice of tie-simplex parameterization for compositional reservoir simulation, in: Proc. of AIChE Annual Meeting, Nashville, TN, 2009. D.Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 15 (1976) 59–64. A. Brøndsted, An Introduction to Convex Polytopes, Springer-Verlag, New York, 1983. G.T. Gilbert, Am. Math. Month. 98 (1) (1991) 44–46. Computer Modelling Group, WinProp Phase Property Program User’s Guide, 2008.