Multiphase Equilibrium Calculation Framework for ... - ACS Publications

Subscriber access provided by UNIV OF BARCELONA

Thermodynamics, Transport, and Fluid Mechanics

Multiphase equilibrium calculation framework for compositional simulation of CO2 injection in low temperature reservoirs Huanquan Pan, Michael Connolly, and Hamdi Tchelepi Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b05229 • Publication Date (Web): 04 Jan 2019 Downloaded from http://pubs.acs.org on January 12, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Multiphase equilibrium calculation framework for compositional simulation of CO2 injection in low temperature reservoirs Huanquan Pan,∗ Michael Connolly, and Hamdi Tchelepi Department of Energy Resources Engineering, Stanford University, Stanford, CA 94305 E-mail: [email protected]

Abstract CO2 injection into an oil reservoir at low temperature can lead to the formation of three hydrocarbon phases. Key challenges persist in the development of phaseequilibrium calculations for three-phase hydrocarbon-CO2 systems. One challenge is the identification of a third equilibrium phase via stability testing of a two-phase mixture. Stability testing requires initial estimates of phase equilibrium ratios (K-values). Existing methods for estimating K-values for the identification of three-phase behavior are unable to detect the whole three-phase region. In addition, many of the initial K-value estimates are redundant; wherein convergence is to the same composition. In this work, we present systematic procedures for single-phase and two-phase stability testing using optimized sets of initial K-value estimates. We show that two-phase stability testing requires to test both of the equilibrium phases in the system. We reveal that three-phase hydrocarbon-CO2 behavior cannot be resolved across the pressurecomposition parameter space with existing methods for the estimates of the initial K-value. We introduce a new initial K-value estimate for two-phase stability testing to

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

resolve three-phase behavior. Our approach significantly improves the reliability and reduces the cost of phase stability testing. Reservoir simulation requires fast and reliable algorithms for stability testing and flash calculations. Failure of the phase equilibrium kernel leads to time-step cuts and the potential cessation of the simulator. We have implemented a combined SS (successive substitution), Newton and trust-region (TR) optimization algorithm for both stability testing and multiphase flash calculations. We have also developed a novel optimization algorithm to solve the multiphase Rachford-Rice (RR) equations. We present the performance of our phase equilibrium framework using nine characterized fluid systems from the literature. Rigorous testing across a wide parameter space demonstrates the robustness of our framework. We show specific instances of quantification of phase equilibrium in challenging scenarios, including near the Stability Test Limit Locus (STLL) in stability testing and the critical region in multiphase flash calculations. The results in this paper address the most prominent difficulties in equation of state (EOS) modeling of CO2 -flooding. In sum, this paper provides the comprehensive detailed algorithms for phase equilibrium calculations for simulation of CO2 -flooding in low temperature reservoirs.

Introduction Injection of CO2 into an oil reservoir at low temperature (less than 322 K) can form up to three hydrocarbon phases: a vapor phase (V), an oil-rich liquid phase (L) and a CO2 -rich liquid phase (L2 ). CO2 has been used for miscible enhanced oil recovery (EOR) in low temperature reservoirs throughout Texas and New Mexico for several decades. 1 In recent years, CO2 -flooding has received additional interest because of the potential to store captured emissions in depleted petroleum reservoirs. 2,3 Miscible displacement with CO2 involves complex interaction of thermodynamics and multiphase flow and requires fully compositional simulation. However, simulation of the CO2 injection process with three hydrocarbon phases

2


Page 2 of 56

Page 3 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


is particularly challenging. To date, multiphase simulation capability does not exist for practical field scale studies involving millions of cells and complex well structures in any commercial simulator. The two-phase representation of three-phase hydrocarbon behavior is standard practice in commercial reservoir simulation. The two approaches are: (i) Only perform two-phase flash; 4,5 (ii) Perform three-phase flash and re-combine the two liquid phases into a single liquid phase. These approaches have limitations. There is no guarantee to remove the entire three-phase region in method (i). 6 Method (ii) may lead to inconsistent flow behavior in simulation. Efficient and robust multiphase equilibrium calculations are essential to enable compositional simulation of CO2 injection at field scale. Information is required on the number of phases present, and the composition and amount of these phases. Although the phase equilibrium calculations for three-hydrocarbon systems have been studied over three decades, two difficulties remain in the development: robust solution of three-phase flash calculations particularly near the critical region, and estimation of equilibrium ratios to initiate phasestability testing. 7,8 We present a framework that systematically addresses both of these challenges.

Numerical Methods Numerical solution of both stability testing and flash calculations usually requires a sequential implementation of algorithms, wherein successive substitution (SS) is paired with a second order method. SS iterations attain a rapid reduction in the magnitude of the residual norm at the beginning of the calculation when the starting point is generally far away from the solution. 7 The SS iterations always advance the result in the direction of the solution. However, SS becomes slow when the iteration is close to the solution. When SS becomes prohibitively slow the Newton Method is typically preferred. Newton’s method generally converges rapidly to the solution with the caveat of a limited radius of convergence. The starting point must be sufficiently close to the solution, or the Newton solver will diverge. 3



Combination of SS and Newton leverages the relative advantage of each method. The combined SS-Newton algorithm first uses SS before switching to Newton’s method. The sequential protocol is generally very efficient, and converges rapidly to the solution in the vast majority of cases. However, the SS-Newton algorithm encounters difficulties when applied to both stability testing and flash calculations. In the case of stability testing, problems occur along the stability test limit locus (STLL), where the tangent plane distance (TPD) is discontinuous, and in some regions above the STLL. 9,10 In the solution of multiphase flash calculations, the numerical difficulty takes place in the near critical region. The SS-Newton method becomes extremely slow when used in these challenging scenarios. If the switch from SS to the Newton method occurs too early, the Newton solver fails and the algorithm must return to advancing the solution using SS iterations. In some cases, the back-and-forth SS-Newton switches occur repeatedly, and a large number of SS iterations are required before Newton’s method is feasible. In our experience, the sole SS-Newton algorithm lacks the robustness and efficiency for phase equilibrium calculations in reservoir simulation. Several authors have proposed globally convergent Newton optimization methods for phase equilibrium calculations. The algorithms used fall into two broad families: (i) Linesearch (LS); (ii) Trust-region (TR). Phase equilibrium implementation of LS was pioneered by Michelsen 7 who used a modified Cholesky decomposition to provide a search direction when the Hessian matrix was not positive definite. 11 Michelsen 8 later recommended the TR method for stability testing and flash calculations for two-phase systems. Trangenstein 12 systematically described the implementation of an LS algorithm for stability testing and flash calculations for two-hydrocarbon-phase systems. Perschke 13 extended Trangensteins algorithm to three-phase hydrocarbon-CO2 systems. Recently, line search with modified Cholesky decomposition has been used for isochoric-isothermal (VT) flash calculations. 14–18 Petitfrere and Nichita 19 implemented a TR method for stability testing and multiphase (up to three) flash calculations. There are many possible trust-region algorithms, 20 a review of which is beyond the scope of this paper. Several researchers 21,22 in the phase behavior

4


Page 4 of 56

Page 5 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


community have implemented Powells dogleg method. 23 Michelsen, 8 Petitfrere and Nichita 19 implemented the More and Sorensen method. 24 However, efficient implementation of these algorithms in stability or flash algorithms is not straightforward. Aside from algorithm complexity, there are questions of efficiency, owing to the overhead of these methods. Indeed, the SS-Newton algorithm is far more effective in most circumstances when applicable for both stability testing and flash calculations. In this paper, we provide a detailed description of a trust-region algorithm in addition to guidance on prudent application. We emphasize that both the LS and TR algorithms are modified Newton methods that converge to a local minimum. These modified Newton solvers are distinct from global minimization algorithms which are briefly discussed below. Derivative-free global minimization algorithms find the global minimum of the underlying function. These algorithms yield the global minimum of the TPD function in stability testing or the minimum Gibbs free energy in the case of multiphase flash calculations. Pan and Firoozabadi 25 used simulated annealing for stability testing and three-phase flash calculations. McKinnon and Mongeau 26 adopted the generic algorithm for calculation of chemical and phase equilibrium. However, derivative-free global optimization methods are prohibitively slow and thus impractical for use in reservoir simulation. Nevertheless, derivative-free methods can be useful to verify the results from the locally convergent algorithms discussed in this research. Quasi-Newton methods have been widely used in application in the phase equilibrium calculations. Quasi-Newton algorithms represent an alternative to full Newton methods, exhibiting super-linear performance in both stability testing and multiphase flash calculations. 9,21,27–30 Nevertheless, for the most difficult flash problems the performance of quasiNewton methods is insufficient, with failure to converge in the critical region after thousands of iterations. 12 The TR or LS approach exhibits greater robustness and efficiency than quasiNewton methods. There are a number of possibilities in the choice of independent variables for flash cal-

5



culations, including conventional and reduced parameters. The primary unknowns are often the component mole numbers, as this allows solution of the phase-split calculation by means of Gibbs free energy minimization. Alternate conventional variables include the logarithm of phase equilibrium ratios (lnK). 29,31 Several variations of reduced variables have also been used in the solution of phase equilibrium calculations. 32–35 Regardless of the choice of independent variables, the underlying challenge of non-convergence persists near the STLL for stability testing and in the critical region for flash calculation.

Solution of Rachford-Rice Equations Phase-split calculations are almost invariably initiated with successive substitution, which requires the solution of the Rachford-Rice (RR) 36 equations at each iteration. The RR equations are solved to obtain phase amounts. This is commonly termed a constant-K-value flash. 37 The number of RR equations is NP − 1, where NP is the number of phases. Solution of the two-phase RR equation is relatively straightforward. However, solving the pair of coupled RR equations in three-phase flash calculations can be challenging owing to the nonlinearity of the objective functions. The inherent nonlinearity can lead to divergence of Newton’s method. Solution of the three-phase RR equations is the subject of several papers in the literature. Leibovici and Neoschil 38 adjusted the step length in each Newton iteration. Michelsen 39 was the first to recast the RR equations into an optimization problem. Leibovici and Nichita 40 later extended Michelsen’s approach by posing the RR equations in a constrained optimization problem with negative flash. Yan and Stenby 41 also extended Michelsen’s approach, introducing a line search into the solution. Taking a different approach, Haugen et al. 31 and Li and Firoozabadi 42 developed a two-dimensional bisection method. Iranshahr et al. 43 used the bisection method for negative RR flash. Okuno et al. 44 further developed the method of Leibovici and Nichita 40 by using a line search algorithm to solve a convex objective function. In addition, Okuno et al. 44 proposed a systematic method to estimate initial phase amounts. The approach involves traversal of the convex polyhedron 6


Page 6 of 56

Page 7 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


formed by the K-value linear constraint equations to determine feasible initial phase amounts for use in multiphase RR equations. A rapid and robust RR solver is an essential element of a multiphase flash algorithm, which almost invariably uses SS when the initial guess is far from the solution.

Strategies for Phase Stability Testing Detection of phase instability requires good initial estimates of phase equilibrium ratios (Kvalues). For a two-phase (V-L) hydrocarbon system, the Wilson correlation 45 and its inverse provide very good initial guesses for stability testing. However, if the system contains more than one liquid phase (L-L2 , or L-L2 -V) it becomes very difficult to make initial estimates of the K-values and there is no guarantee of detecting phase instability. Michelsen 46 first suggested using a set of initial K-value estimates for two-phase stability testing, which included: (i) Near-pure corner point trial phase compositions corresponding to each component in the mixture; (ii) an ideal state in the trial phase; (iii) an average composition of two-phase compositions in the trial phase. Stability testing in three-phase systems used the suggested initial K-value estimates for several years. 13,28 Li and Firoozabadi 29 extended Michelsen’s approach 46 with a set of initial K-value estimates comprising the Wilson correlation and its inverse, the cubic root of the Wilson correlation and its inverse, and corner point estimates corresponding to 90%-mole fraction of each component in the trial phase. We show that even when this large set of proposed estimates is used, failure to detect two-phase instability occurs for several fluids. Li and Firoozabadi 29 used a total of Nc + 4 initial K-value estimates in both single-phase and two-phase stability testing for a wide range of fluids, including low temperature sour gas mixtures containing H2 S and asphaltene-containing mixtures. For hydrocarbon-CO2 mixtures, we show that most of the Nc + 4 initial estimates are redundant, and that only a few are required. In this research, we propose a smaller set of K-value estimates. An overlooked issue in phase stability testing is test phase selection when more than 7



one phase is present. Detection of a third equilibrium phase requires stability testing of a two-phase system. Perschke 13 and Li and Firoozabadi 29 noted that test phase selection does not influence stability testing, the outcome of which will be the same based on Gibbs free energy analysis of the equilibrium phases by Baker et al. 47 However, the tangent plane criterion developed by Baker et al. 47 requires the global minimum of the TPD function. Since stability test calculations are solved using local minimization methods, it is not safe to assume that stability testing is insensitive to the selection of the test phase. We show that given the same initial K-value estimate for a two-phase state, stability testing can yield different results depending on the phase tested for instability. For instance, a two-phase mixture may appear stable when testing the lighter phase, and unstable when testing the heavier phase. The implication is that detection of three-phase equilibrium is unreliable when only one of the phases in a two-phase mixture is subject to stability testing.

Implementation and Overview The Stanford ADGPRS (Automatic Differentiation based General Purpose Research Simulator) provides a flexible framework for the development of multiphase compositional simulation capability. 48 The purpose of this study is to develop rapid and reliable multiphase equilibrium calculations for three-phase hydrocarbon-CO2 systems. Our objective is to provide the phase equilibrium kernel for ADGPRS, and simulate three-phase flow in CO2 -floods of low temperature reservoirs. In this paper, we focus on developing a comprehensive framework for multiphase equilibrium calculations in hydrocarbon-CO2 systems. We develop a sequential SSI-Newton-TR numerical iteration schema to achieve high efficiency and robust convergence in stability testing and multiphase flash calculations. We also develop a novel RR solver by using a TR optimization method. We use the minimum number of initial K-value estimates for reliable single-phase stability testing. Previously proposed initial Kvalue estimates for two-phase stability testing are inadequate for identification of three-phase equilibrium. We introduce a new initial K-value estimate to supplement existing estimates 8


Page 8 of 56

Page 9 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in two-phase stability testing. We also optimize the number of initial K-value estimates for two-phase stability testing. We tested our solution framework extensively with nine characterized fluids taken from the literature. Most of the characterized fluids are from reservoirs in Texas and New Mexico. For each mixture, we constructed the entire phase envelope by performing phase equilibrium computations over a wide range of pressures and injected CO2 mole fractions. The results show that our protocol for phase equilibrium computations correctly resolves all single-phase, two-phase and three-phase regions. The paper proceeds as follows. First, we formulate the phase stability testing problem, as well as the two-phase and three-phase flash calculations. Next, we describe the numerical algorithms to solve the stability testing and flash calculations. We then propose a new method to estimate the initial K-values used in phase stability testing. We present and discuss the phase diagrams and the numerical results for tested fluids. Finally, we summarize and draw conclusions.

Formulation Single-Phase Stability Testing Given a mixture with composition zi (i = 1, .., Nc ) at temperature T and pressure P , singlephase stability testing determines whether the mixture is unstable as a single phase. The stability testing calculation is based upon minimization of the tangent plane distance (TPD) function: TPD = 1 +

Nc X

Yi

V L ˆ ˆ lnφi (Yi ) + lnYi − lnφi (zi ) − lnzi − 1

(1)

i=1

where Yi is the mole number of component i in the trial phase, 46 and φî is the fugacity coefficient which is evaluated using the Peng-Robinson equation of state (PR-EOS). 49 If the optimal TPD is negative, then the fluid is unstable. Otherwise it is stable. If the fluid is

9



Page 10 of 56

unstable, the phase equilibrium ratios, Ki = Yi /zi (i = 1, 2, ..Nc ) from the stability test are used as the initial guess for the ensuing two-phase flash calculation.

Two-Phase Stability Testing Given a two-phase mixture with equilibrium phase compositions of xi (heavy phase, L) and yi (light phase, V) (i = 1, ..Nc ) at temperature T and pressure P , two-phase stability testing checks for instability of the two-phase system. A three-phase flash calculation is required upon detection of two-phase instability. A two-phase stability test calculation is identical to those for single-phase stability testing, except in Eq. 1 the feed composition zi is replaced with the composition of one of the two equilibrium phases (xi or yi ).

Two-Phase Flash Calculation Two-phase (V and L) flash calculations are performed through minimization of the Gibbs free energy at given feed composition zi (i = 1, ..Nc ), pressure P , and temperature T ,

¯ = G/RT = G

Nc X

nVi lnfîV (nVi )

i=1

+

Nc X

nLi lnfîL (nLi )

(2)

i=1

where fî refers to the fugacity of component i. The optimization is subject to the following material balance constraints:

nLi + nVi = zi

(i = 1, .., Nc ).

(3)

The independent variables are either component mole numbers of light phase (V ) nVi or the heavy phase (L) nLi (i = 1, .., Nc ).

10


Page 11 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Three-Phase Flash Calculation Minimization of Gibbs free energy is the basis for our three-phase flash implementation, as for the two-phase flash. For a three-phase system the Gibbs free energy at given feed composition zi (i = 1, ..Nc ), pressure P , and temperature T is,

¯ = G/RT = G

Nc X

nVi lnfîV (nVi )

i=1

+

Nc X

nLi lnfîL (nLi )

+

i=1

Nc X

nLi 2 lnfîL2 (nLi 2 ).

(4)

i=1

The optimization is subject to the following material balance constraints:

nLi + nVi + nLi 2 = zi

(i = 1, .., Nc ).

(5)

We use the heavy liquid phase L as the reference phase. The independent variables are the component mole numbers of the vapor phase nVi and the light liquid phase nLi 2 (i = 1, .., Nc ).

Numerical algorithms Numerical Algorithms for Stability Testing and Flash Calculations Stability test calculations proceed via minimization of the TPD function in Eq. 1. For ¯ in Eq. 2, subject to the material balance the two-phase flash calculation, we minimize G ¯ in constraints (Eq. 3). Similarly, the three-phase flash calculation involves minimization of G Eq. 4 with the material balance constraints in Eq. 5. We use the widely accepted numerical protocol in phase equilibrium calculations; Newton’s method is employed to attain rapid convergence after first using successive substitution. However, we add robustness in the form of a trust region solver, which we have implemented to provide convergence in the case of difficult problems for both stability testing and flash calculations. For optimum efficiency, our implementation leverages the relative strength of the SS, Newton and TR algorithms. Namely, solution of the stability and flash problems proceeds in a sequential format: (i) SS

11



Page 12 of 56

steps provide initial convergence far from the solution; (ii) the Newton method yields rapid convergence once the iteration point is close to the solution; (iii) the TR solver is used when a difficult problem is encountered and convergence with Newton is impractical. In summary, we use a sequential SS-Newton-TR iteration algorithm for minimization of Eqs. 1, 2 and 4 (see Fig. 1 for the algorithm flowchart). The major application of the TR solver is restricted to cases for which the standard SS-Newton solution protocol is inadequate, such as where the TPD function is close to the ¯ function is in the critical region (flash). In these scenarios, STLL (stability testing) and the G Newton’s method may encounter infeasible conditions. The infeasible condition in Fig. 1 within a Newton iteration is

Yi < 0

or

T P Dnew > T P D,

(6)

for stability testing using Eq. 1, and

nVi < 0

or

nLi < 0

or

¯ new > G, ¯ G

(7)

for two-phase flash calculations using Eq. 2, and

nVi < 0

or

nLi < 0

nLi 2 < 0

or

or

¯ new > G ¯ G

(8)

for three-phase flash calculations using Eq. 4. For stability testing, assessment of convergence is with reference to the residual norm of the gradients of Eq. 1: v u Nc 2 uX V L ln Yi + ln φî (Yi ) − ln zi − ln φî (zi ) . =t

(9)

i=1

For two-phase flash calculations, the residual norm of the fugacity equations is the measure

12


Page 13 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of convergence: v u Nc 2 uX =t fîV − fîL .

(10)

i=1

Similarly, for three-phase flash calculations the norm of interest is, v u Nc 2 2 uX L2 t V L L ˆ ˆ ˆ ˆ fi − fi + fi − fi . =

(11)

i=1

In our extensive testing, we found that the computational performance is insensitive to the switching criterion switch value from SS to Newton in both stability testing and flash calculations. For stability testing we switch from SS to Newton once switch =0.1. For twophase and three-phase flash calculations, we use a tighter threshold of switch =0.01 before switching to Newton’s method from SS. The selection of the switch values is discussed in more detail later. We also place a limit on the maximum SS iterations permitted. An upper limit of 20 iterations prevents excessive SS iterations, which can occur close to the STLL in the stability testing and the critical point in flash calculations. When Newton’s method encounters an infeasible condition, we discard the results from the Newton iterations and initiate the TR solver using the K-values from the last SS iteration (see Figure 1). Detailed descriptions of SS and the Newton method as applied to phase equilibrium calculations are available in the literature. See Michelsen and Mollerup 50 for details. A brief description of our implementation of a TR optimization algorithm is provided herein. At each iteration, a TR algorithm comprises two major steps: (i) Update the TR radius ∆K ; (ii) compute the step-size sK for the restricted step within ∆K . The second operation of the algorithm is called the trust-region sub-problem. Several methods have been proposed for the solution of the sub-problem, 51 a comprehensive discussion of which are beyond the scope of this paper. We restrict our review to trust region implementations in phase equilibrium calculations.

13



Page 14 of 56

Two algorithms have been used to address the TR sub-problem for phase-equilibrium computations. The first is the dogleg method. 21–23 The second is the iterative solution method originally proposed by More and Sorenson, 24 which has been implemented by Michelsen 8 and more recently by Petitfrere and Nichita. 19 We have implemented the iterative method for the solution of the TR sub-problem, as we believe it to be more suitable for stability testing and flash calculations. We closely follow the algorithmic description provided in Nocedal and Wright, 51 Figure 2 provides a schematic description of the algorithm used to compute ∆K . In Fig. 2, ρk is defined as

ρk =

f (xk + sk ) − f (xk ) , m (xk + sk ) − m (xk )

(12)

where f (xk ) is the objective function. The model function m (xk ), is defined as 1 m (xk ) = f (xk ) + gk sk + sTk Hsk , 2

(13)

where g is the gradient vector and H is the Hessian matrix of the objective function f (xk ). Note that if a decrease in the objective function is not realized (i.e. ρk < 0) we use a single SS iteration to compute xk+1 from xk . In contrast, the standard TR algorithm does not update the independent variables (i.e. xk+1 = xk ) when ρk < 0 (see Figure 2). Petitfrere and Nichita 19 first suggested the SS update to advance the solution in the case of an increase in the objective function. We have found the SS update greatly improves the performance of the TR solver in difficult situations. The TR sub-problem requires the solution of

(H + λI) s = −g,

(14)

subject to the constraint ksk2 ≤ ∆ to obtain the step s. Here, λ is a scalar and λ ≥ 0. Fig. 3 provides an algorithmic description of the solution to

14


Page 15 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the TR sub-problem. Note that if the Hessian matrix H is not positive definite, it is necessary to calculate its smallest eigenvalue λ1 for the purposes of re-setting λ so that the matrix (H+λI) is positive definite. Calculation of the minimum eigenvalue is an expensive operation. However, the TR-algorithm is always preceded by SS iterations, and H is generally positive definite. For this reason, it is rarely necessary to calculate λ1 . The iterative procedure to calculate s generally converges within two or three iterations.

Solution of the Three-phase Rachford-Rice (RR) Equations At each SS iteration in a three-phase flash calculation, it is necessary to solve the following RR equations: Nc X i=1

Nc X i=1

zi (KiV − 1) = 0, 1 + V (KiV − 1) + L2 (KiL2 − 1)

(15)

zi (KiL2 − 1) =0 1 + V (KiV − 1) + L2 (KiL2 − 1)

(16)

which yield the vapor phase amount V and light liquid phase amount L2 . In Eqs. 15 and 16, the superscripts V and L2 represent the vapor and light liquid phase respectively. Solution of the inherently nonlinear RR equations is difficult, particularly in the case of a poor initial estimate for the phase amounts. Michelsen 39 was the first to pose solution of the RR equations as a minimization problem. Leibovici and Nichita, 40 Okuno et al., 44 and Yan and Stenby 41 further extended the Michelsen approach to negative flash. We use the objective function, proposed by Okuno et al., 44 as follows:

FRR (V, L2 ) =

Nc X

−zi ln 1 + V (KiV − 1) + L2 (KiL2 − 1) .

i=1

15


(17)


Page 16 of 56

Compared with our approach, Leibovici and Nichita 40 use an objective function with an extra term: summation of all phase amounts in Eq. 17. We use a TR method for optimization of phase amounts in Eq. 17. The constraints in the optimization are: 40

1 + V (KiV − 1) + L2 (KiL2 − 1) > 0 (i = 1...Nc ).

(18)

The phase amounts V and L2 are permitted to become negative during the iterative solution of Eq. 17, provided the constraints of Eq. 18 hold. 40 The residual is the Euclidean norm of Eqs. 15 and 16. In our implementation, we adjust the relaxation parameter t so that the constraints t 1 + V (KiV − 1) + L2 (KiL2 − 1) > 0 (i = 1...Nc ) are satisfied upon initialization and within each iteration. With appropriate adjustment of t, our TR algorithm is both robust and insensitive to the initial estimate for phase amounts. Details regarding derivation and estimation of the relaxation parameter t are provided in a separate publication. 52 In this research, 52 we provide a framework to calculate multiphase (up to four phases in examples) Richford-Rice equations using TR-based optimization algorithms in multiphase compositional and thermal simulation.

Initial Estimates of Phase Equilibrium Ratios (K-Values) for Stability Testing Conventional approaches to phase stability testing use several initial K-value estimates, derived from the Wilson correlation, ideal gas law and corner point estimates corresponding to near pure trial phases: 7,29

KiW ilson , 1/KiW ilson ,

q 3

KiW ilson , 1/

q 3

KiW ilson , KiP ure−j (j = 1...Nc ), KiIdeal , KiAverage .

16


(19)

Page 17 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The correlations of these estimates are

KiW ilson = Pci exp[5.37(1 + ωi )(1 − Tci /T )]/P, KjP ure−j = 0.9/zj , KiP ure−j = 0.1/(Nc − 1)/zi

(i 6= j, i = 1, ...Nc ),

(20)

(21)

KiIdeal = φî (zi ),

(22)

KiAverage = (xi + yi )/2.

(23)

Note that KiP ure−j (j = 1, ..Nc ) consist of Nc estimates. The KiAverage estimate is only required for two-phase stability testing. Hence, xi and yi refer to the equilibrium phase compositions of the two-phase mixture. In total, Nc + 5 estimates were proposed for singlephase stability testing and Nc + 6 for two-phase stability testing. In the case of strictly two-phase vapor-liquid systems, KiW ilson and 1/KiW ilson provide excellent initial guesses for stability testing. However, complexity arises for the estimation of initial K-values when the system contains two-liquid phases (either liquid-liquid or vapor-liquid-liquid). For mixtures comprising more than one liquid phase there is no universal criteria for the estimation of K-values. Li and Firoozabadi 29 proposed using the following Nc + 4 estimates for both single-phase and two-phase stability testing:

{KiW ilson , 1/KiW ilson ,

q 3

KiW ilson , 1/

q 3

KiW ilson , KiP ure−j (j = 1...Nc )}.

(24)

We will show that these estimates are both incomplete (i.e. fail to detect two-phase instability at times) and redundant for the CO2 -hydrocarbon fluids in this research. The redundancy is manifest in the identical trial phase compositions computed using stability testing with different initial estimates. In our experience, for a three-phase system the initial estimates should depend on the specific type of fluid. For example, initial K-value estimates for a three-phase vapor-liquid-water system should be different from the initial estimates for three-

17



Page 18 of 56

phase hydrocarbon-CO2 systems. Since water is essentially immiscible with hydrocarbons, the distinctive behavior of hydrocarbon-aqueous-steam systems yields a priori estimates of initial K-values. 53 Hydrocarbon-water systems require far fewer initial K-values for stability testing than highly miscible hydrocarbon-CO2 systems.

Initial K Estimates for Single-phase Stability Testing For single-phase stability testing, we use only five initial K-value estimates:

{KiW ilson , 1/KiW ilson ,

q 3

KiW ilson , 1/

q 3

KiW ilson , KiCO2 }.

(25)

We have found that adding extra K-value estimates to those in Eq. 25 does not yield new trial phase compositions in the stability test calculations. We calculate TPD values for all five initial estimates in Eq. 25. The converged K-values yielding the lowest negative TPD serves as the initial guess for the ensuing two-phase flash calculation.

Initial K Estimates for Two-phase Stability Testing Stability testing of a two-phase fluid with all Nc + 6 published initial estimates results in occasional failure to detect three-phase behavior. Detecting instability of the two-phase mixture is unreliable, even when both the L and V phases are tested. New initial K-value estimates are required for consistent identification of two-phase instability. Based on the phase behavior of three-phase hydrocarbon-CO2 systems, we propose a supplementary initial pure CO2 estimate: CO2 −95 KCO = 0.95/zCO2 , KiCO2 −95 = 0.05/(Nc − 1)/zi 2

(i 6= CO2 , i = 1, ...Nc ).

(26)

In Eq. 26, the superscript 95 refers to the 95% mole fraction of CO2 in the trial phase (note that the CO2 concentration is 90% in the trial phase in the initial estimate KiCO2 . We are

18


Page 19 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


able to identify the existence of previously un-detected three-phase points in all our test cases using the new estimate. To detect instability of a two-phase mixture consistently, it is necessary to test both the heavy L-phase and the light V-phase. For stability testing of the L-phase, we propose the following set of initial estimates:

{KiCO2 , KiW ilson ,

q 3

KiW ilson , KiCO2 −95 }.

(27)

Meanwhile, we use a different set of initial estimates for stability testing of the V-phase:

{KiCO2 , 1/

q 3

KiW ilson , KiCO2 −95 }.

(28)

In total, we require seven stability test calculations to test for the appearance of a third phase. The detection of a negative tangent plane distance indicates the instability of a phase. For a strictly two-phase vapor-liquid system two-phase flash is performed upon the first detection of instability. Once a stability test yields a negative TPD, the two-phase flash calculation proceeds using the K values from the stability calculation as the initial estimate of phase compositions. To be clear, it is not necessary to perform additional stability test calculations using other initial K-value estimates. However, for three-phase flash calculations, the result of a two-phase stability test corresponding to a negative TPD is not necessarily a good starting point for the flash calculation. In our experience, this is a key difference between stability testing in two-phase and three-phase systems. We attribute this difference to the complexity of three-phase mixtures. For a system in which a maximum of two-phases may exist, it is not always necessary to perform stability tests using both KiW ilson and 1/KiW ilson as initial estimates. Once a stability test yields a negative TPD, we immediately perform a two-phase flash calculation using the K values from the stability calculation as an initial estimate for the phase split calculation. The approach works very well for two-phase vapor-liquid systems. 19



Theoretically, the same approach might be implemented for a mixture exhibiting three-phase behavior; cease two-phase stability testing upon first detection of a negative TPD, and use the K values from the stability test to initiate the three-phase flash calculation. Due to the complexity of three-phase systems, on rare occasions non-trivial and very small negative TPD values could occur in two-phase stability testing, attributable to convergence to a local minimum. In this scenario, the K values from the stability testing calculation constitute an infeasible set, and lead to failure of the three-phase flash calculation. Our approach is to proceed with two-phase stability testing until all seven initial K-value estimates listed in Eqs. 27-28 are used. If each of the calculated TPD values is larger than zero, the fluid is in a two-phase state. Otherwise, we select the K values associated with the lowest TPD value and use them as the initial guess for the ensuing three-phase flash calculation. Occasionally, the three-phase flash calculation fails. When this occurs, we return to the single-phase stability test results, select the set of K values corresponding to the next lowest TPD, and perform a new two-phase flash calculation. We then proceed with two-phase stability testing and a three-phase flash. Figure 4 plots the flow-chart diagram for the algorithm of the multiphase equilibrium calculations. Note that the failed three-phase flash calculation is generally from the two-phase flash result from which the Gibbs free energy is not lowest. We further discuss this issue in the following section.

Results and Discussion In this section, we discuss the results of extensive testing of our framework for phase equilibrium calculations involving CO2 . We tested nine fluids published from the literature. Most of these fluids are from CO2 EOR projects in West Texas and New Mexico. Tables A1-A9 describe the thermodynamic properties of each fluid. The non-zero binary interaction coefficient kij values are also listed in these tables. Table A10 lists the reservoir temperatures of each of the mixtures. These data are directly obtained from the publications of Khan et

20


Page 20 of 56

Page 21 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


al. 54 (BSB, JEMA, NWE and OilG), Li and Firoozabadi 29 (Maljamar Reservoir Oil (MRO), Maljamar Separator Oil (MSO) and OilB) and Pan et al. 5 (ResA and ResB).

P-x Phase Diagrams We developed a comprehensive test suite using Matlab. Using the test suite, we generate input files with the test resolution specified by the user. The test suite also has options for post-processing results. For each fluid, we divide the P-x (pressure vs. injected mole fraction of CO2 in the mixture) diagram into three regions with different resolutions: (i) a three-phase region across which we perform phase equilibrium calculations at high resolution; (ii) a twophase liquid-liquid (LL2 ) region; (iii) a single-phase region where the concentration of CO2 is very high. The coordinates of each tested point are xCO2 (injected mole fraction) on the X-axis and pressure (bar) on the Y-axis. Figs. 5a-5i show the P-x phase diagrams for the nine fluids in our test suite. Pressures and compositions were selected so that the three-phase region is visible for each fluid system. We consistently navigate phase boundaries. Note that we do not distinguish between vapor and liquid phase states using phase labels, nor do we classify the two-phase regions (VL, VL2 or LL2 ). Phase labeling is an outstanding issue in the study of fluid phase behavior of hydrocarbon-CO2 system at low temperatures, 13 and is not discussed in this research. Additional P-x phase diagrams are presented to resolve the liquid-liquid-phase region for all tested fluids. Figures 6a-6c display P-x phase-diagrams for NWE, OilG and ResA fluids respectively. Of particular interest is the presence of a singlephase CO2 -rich liquid (L2 ) in the upper right corner of each of these phase diagrams. Our test results show that the CO2 -rich single phase exists at very high CO2 concentrations for all of the fluids tested. Figs. 7a-7c show the P-x diagram at the very high CO2 concentrations for the NWE, OilG, and MSO fluids, respectively. MSO is a separator crude oil, rather than a reservoir fluid. However, at elevated pressures CO2 still dissolves MSO to form a single-phase (See Figure 7c). Although the single-phase CO2 -rich regions are very small and difficult to see in a P-x diagram, the identification of single-phase behavior is important 21



in the context of compositional reservoir simulation. In miscible gas floods, injection of high concentrations of CO2 may lead to single-phase behavior near injection wells. Li and Firoozabadi 29 only reported the detection of a CO2 -rich single phase for BSB. Our phase equilibrium calculation framework is able to detect the CO2 -rich liquid phase in all of the fluids studied in this research. For verification purposes, we also performed phase equilibrium calculations and constructed P-x phase diagrams using all of the initial K-value estimates proposed in the literature. For single-phase stability testing there are a total of Nc + 4 estimates from the literap p ture: KiW ilson , 1/KiW ilson , 3 KiW ilson , 1/ 3 KiW ilson , KiP ure−j (j = 1...Nc ). For two-phase stabilp p ity testing there are a total of Nc + 7 estimates: KiW ilson , 1/KiW ilson , 3 KiW ilson , 1/ 3 KiW ilson , KiP ure−j (j = 1...Nc ), KiCO2 −95 , KiIdeal , KiAverage . We test both the V and L phases, resulting in 2 x (Nc + 7) trials for each two-phase stability test. We refer to stability testing performed using the full spectrum of proposed initial K-values as the reference method. Our implementation uses a smaller set of K-values. Herein, we refer to our approach as the optimized method. When using our optimized method we require five calculations for single-phase stability testing and seven calculations for two-phase stability testing. For each of the nine fluids in this study, we tested the reference method and the optimized method using the same input files. The implementations are identical, aside from the initial K-values used in stability testing. We compared the phase diagrams generated using each of the two approaches to verify the results from the small set of initial K-values estimates in the optimized method.

Initial K-value Estimates for Stability Testing In the case of single-phase stability testing, computation of redundant trial phase compositions occurs when using corner point initial K-value estimates other than those derived from the injection component ( CO2 ). Redundant trial phase compositions are those that are repeatedly attained when initiating the stability test with different sets of initial K-value estimates. In our experience, non-CO2 corner-point initial K-value estimates lead to trial 22


Page 22 of 56

Page 23 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


phase compositions identical to those already computed using either KiW ilson , 1/KiW ilson or KiCO2 . For this reason, we exclude non-CO2 corner point initial K-value estimates from single phase stability testing. Occasionally, the only means to detect instability of the mixture is by using initial K-value estimates derived from the cubic root of the Wilson correlation and p p p its inverse 3 KiW ilson , 1/ 3 KiW ilson . For example, 3 KiW ilson is only initial estimate to detect mixture instability of the MSO fluid at points of (0.682 70.22) and (0.682 70.36). Through extensive testing, we have found that five sets of initial K-value estimates work very well. Perschke, 13 Li and Firoozabadi 29 suggested that selection of the test phase (V or L) does not affect the two-phase stability test result. Their explanation was that the TPD function in Eq. 1 is the same value, regardless of the phase selected for testing, due to the equality of lnzi − lnφî (zi ), i.e. lnxi − lnφî (xi ) = lnyi − lnφî (yi ) at equilibrium in the two-phase flash. However, complication arises in the calculation of mole numbers Yi = Ki zi (zi = xi for test L phase, and zi = yi for test V phase) in Eq. 1 and (xi 6= yi ). With the same initial K-value estimate, the initial Yi value differs, and depends on the composition of the test phase. Since we use a locally convergent method to minimize the TPD function in Eq. 1, the different initial Yi value generated using the L or V phase can lead to different trial phase compositions in two-phase stability testing. If we use a global optimization method for minimization of the TPD function Eq. 1, the result of two-phase stability testing will be insensitive to the initial Yi value, and the assumption made by Perschke 13 and Li and Firoozabadi 29 thus holds. Unfortunately, global minimization algorithms are impractically slow for use in reservoir simulation. Given our desire to use a locally convergent algorithm for stability testing, it is imperative we consider the influence that the test phase selection has on the results of two-phase stability testing. Li and Firoozabadi 29 suggested a total of Nc + 4 initial estimates KiW ilson , 1/KiW ilson , p p 3 KiW ilson , 1/ 3 KiW ilson , KiP ure−j (j = 1...Nc ) for both single-phase and two-phase stability testing. We use this set of initial K-value estimates to test heavy L phase, light V phase and both L and V phases in two-phase stability testing. Note that we need 2x(Nc + 4)

23



stability test calculations if both V and L phases are tested for a two-phase stability testing. Figures 8-10 illustrate the results of detecting three-phase region for BSB, JEMA and NWE fluids. These figures show that only testing one phase (L or V) leads to an incomplete three-phase region in the P-x diagram. The two-phase stability testing yields better results when testing the light V-phase for instability, rather than the L-phase. Use of the V-phase results in fewer missing three-phase points for the three fluids tested here. This contradicts the recommendation of Li and Firoozabadi 29 who suggested that only the L phase should be tested. We also observed that the initial estimates KiIdeal and KiAverage are not useful for stability testing of the three-phase hydrocarbon-CO2 systems in this study. The value of the TPD function is always non-negative from these estimates. Li and Firoozabadi 29 also excluded them. Consistent identification of all three phase regions is not possible when using the Nc + 4 p p initial estimates KiW ilson , 1/KiW ilson , 3 KiW ilson , 1/ 3 KiW ilson , KiP ure−j (j = 1...Nc ) for phase stability testing, even if both the V and L phases are tested. Among the nine fluids tested in this study, stability testing with the aforementioned Nc + 4 estimates fails to identify three-phase behavior for the JEMA (Figure 9c), NWE (Figure 10c), MRO and ResB fluids. Figure 11a shows inconsistent identification of three-phase behavior for the MRO fluid. These points are present at low concentrations of injected CO2 . For the ResB fluid, failure to detect three-phase behavior occurs at high concentrations of injected CO2 (Figure 11b). Table 1 lists some of the pressures and compositions at which the method fails to identify three-phase behavior for these four fluids. We detect instability of the two-phase state for each of the points in Table 1 using our proposed new initial estimate KiCO2 −95 , with either the L or V phase taken as the test phase (see Figs. 5b, 5c, 5e and 5i). With our proposed set of initial K-value estimates, we correctly determine all phase states for the nine fluids. We expect that our stability testing strategy will lead to consistent identification of the phase state in hydrocarbon-CO2 systems of interest. However, we cannot make this guarantee, given that estimation of K-values to initiate stability testing is

24


Page 24 of 56

Page 25 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


inherently heuristic. In our experience, varying the mole-fraction of CO2 in the trial phase yields the biggest improvement in the ability to detect the appearance of a third phase. In the event that the initial estimates proposed in this paper are unable to resolve regions of three-phase behavior, we recommend an additional set of K-values that modifies the CO2 mole fraction of the trial phase.

Solution Strategies for Three-phase Systems One of the challenges encountered in three-phase mixtures is the occurrence of several negative TPD values during single-phase stability testing. The conventional approach is to perform two-phase flash with the set of K-values that corresponds to the lowest negative TPD. However, the two-phase flash calculation initiated with this set of K-values does not ¯ Occasionally, a lower G ¯ results from a twoalways yield the lowest Gibbs free energy, G. phase flash initiated with K-values corresponding to the next lowest negative TPD value. There are implications for the ensuing two-phase stability analysis. First, convergence of the ¯ results in an unstable two-phase state. A negative two-phase flash to a local minimum G TPD is then realized in two-phase stability testing. Finally, the three-phase flash calculation fails. In this circumstance, the solution is to perform an additional two-phase flash calculation using the K-values corresponding to the next lowest negative TPD, as generated from single-phase stability testing. Upon convergence of the two-phase flash we test the two-phases for stability (see Figure 4). In the majority of cases, the system is two-phase. Li and Firoozabadi 29 made the same observation. However, on occasion, the system is threephase. Table 2 provides examples of the difficult-to-resolve three-phase behavior. A relevant example is the point x = 0.793, P = 69.46 bar for the NWE fluid. Table 3 lists the results for single-phase and two-phase stability testing, and for two-phase and three-phase flash. Following stability testing, we perform a two-phase flash with the K-values corresponding to the minimum negative TPD. The two-phase flash result is unstable, and the ensuing three-phase flash fails. We perform a second two-phase flash calculation using the K-values associated 25



with the second lowest negative TPD. The resulting two-phase state possesses a lower Gibbs free energy than that obtained with the first two-phase flash. Two-phase stability testing indicates the presence of a third equilibrium phase. Finally, we perform a three-phase flash to quantify the equilibrium state of the mixture. An alternative approach for rigorous quantification of three-phase equilibrium is to perform multiple two-phase flash calculations. The rationale here is that we can avoid failed three-phase flash calculations via stability testing of the two-phase state with the minimum Gibbs free energy. We tested this alternate approach by simultaneously performing two-phase split calculations with the K-values corresponding to the two smallest TPD values obtained in single-phase stability testing. We refer to this approach as double-2P-flash. In the vast majority of cases the flash results attained from each of the two-phase flash calculations are identical. In the event that the two-phase flash results differ, the two-phase mixture with ¯ is selected and two-phase stability testing then proceeds. We found that this the lower G approach greatly reduced the incidence of false two-phase instability and failed three-phase flash calculations. We do not recommend the double-2P-flash approach for two reasons. First, the computational cost of a failed three-phase flash calculation is typically quite low. The reason is that failure of the three-phase flash generally occurs within the first few SS iterations, owing to a negative phase amount computed via solution of the Rachford-Rice equations. Second, the computational cost of performing additional two-phase flash calculations is very high in the double-2P-flash approach; the method requires an extra two-phase flash calculation at every point located in the two-phase or three-phase region when single-phase stability testing yields more than one negative TPD. For example, consider the two-phase diagram for the NWE fluid in Figure 6a. In total, there are 93,600 test points. Execution time is 35 seconds for the standard method (Figure 4), and 61 seconds for the double-2P-flash approach. In the calculation of the three-phase diagram (Figure 5c) which consists of 630,450 test points, the execution time is 413 seconds for the standard approach and 420 seconds for the double-

26


Page 26 of 56

Page 27 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2P-flash approach. In the three-phase region, the majority of the computational time is attributable to the three-phase flash calculation itself. Therefore, in the three-phase region an additional two-phase flash calculation is not a significant overhead. However, the cost of an additional flash in the two-phase region is very high. The three-phase region is much smaller than the two-phase region in the P-x diagram. In reservoir simulation, three-phase behavior is present in only a very small fraction of the grid blocks at a given time step. The double-2P-flash method has a strong theoretical rationale: Identification of threephase equilibrium via stability testing of the two-phase state with the lowest Gibbs free energy. Practically, the approach costs more time and provides little improvement in twophase stability testing and three-phase flash calculations. Therefore, we recommend the standard approach (Figure 4) instead of the double-2P-flash method.

Numerical Solution of Stability Testing and Flash Calculations Our implementation of the numerical algorithm in Figure 1 is both robust and efficient. Through rigorous testing, we covered tens of millions of sample points for the nine fluids in this study. We have not encountered a single point of failure. A given test point requires five calculations for single-phase stability testing, corresponding to the five initial equilibrium ratio estimates. Upon detection of single-phase instability, we require an additional seven calculations for two-phase stability testing. Stability testing calculations consume most of the computational time. An efficient and robust algorithm for stability testing forms the basis of high performance phase-equilibrium computations. We improved the performance √ of our stability testing algorithm by using the iterative variables αi = 2 Yi (i = 1, .., N c) proposed by Michelsen, 46 in both the Newton and TR solvers. In addition, we implement a check for early detection of convergence to the trivial solution. 46 The trivial solution is the solution for which the composition of the trial phase converges to that of the test phase and Ki = 1 (i = 1, .., N c). We terminate the iterations upon detection of convergence to the trivial solution. This greatly reduces total computational time, as convergence to the 27



trivial solution is very common. We strongly advocate early detection and termination of the trivial solution in any stability-testing algorithm. Convergence of stability test calculations is difficult near the STLL. 9,10 In our testing, we frequently encountered points close to the STLL, particularly in two-phase stability tests. We performed numerical experiments as follows. First, we fixed the maximum number of SS iterations to 2000. Generally the residual was reduced to a tight tolerance ( < 10−4 ) after 2000 SS iterations. At this point, we switched to Newton’s method. Despite a tight switching criterion, the Newton method failed to converge for many points. In another numerical experiment, we set switch < 10−5 as our criterion. Again, we found that Newton’s method still failed to converge for some test points. Figure 12 shows the results of two-phase stability test calculations at point (0.603, 82.4) for OilG. The initial K-value estimate is KiCO2 and the light V-phase is the test phase. Figure 12a shows the iteration results using SS-Newton with a switch criterion of switch < 10−5 . The Newton method is initiated after 3304 SS iterations. The residual oscillates and fails to converge after 15 Newton iterations (the maximum permissible Newton iterations). We require an additional 5425 SS iterations to converge to a tolerance of < 10−6 , amounting to a total of 8729 SS iterations for convergence. In contrast, we are able to achieve rapid convergence using our SS-Newton-TR algorithm (see Figure 1) with a default switch criterion of switch < 0.1. Figure 12b shows the results of the stability testing calculation using our numerical implementation. Just two SS iterations are adequate to reach the switch threshold, at which point we switch to Newton’s method. After 3 iterations Newton’s method encounters an infeasible condition, and we switch to our TR solver. Convergence is attained after just 11 TR iterations, including one SS update (see Figure 2 for the step update in TR algorithm). Note that we use the K-values from the last SS iteration to initiate the TR solver. We observed that our TR solver converges in a maximum of 15 iterations for any difficult point close to the STLL. The TR solver displayed this consistency for stability test calculations involving each of the nine fluids in our test suite. Our approach is an improvement on

28


Page 28 of 56

Page 29 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that of Li and Firoozabadi, 29 who assumed that the solution to the stability test is trivial if the iteration number exceeds 1000 using the BFGS-quasi-Newton method. In reservoir simulation, we encounter points close to the STLL frequently and cannot afford to perform 1000 quasi-Newton iterations. For two-phase flash calculations, we observed that a combination of SS and Newton iterations is inadequate for convergence at some points. For example, consider the point (0.747, 77.06) for OilG. At this point, the system is three-phase. We performed the two-phase flash calculation using the initial K values obtained from the single-phase stability testing calculation. For illustrative purposes, we set switch = 10−5 as our criterion for switching from SS to Newton. After 1010 SS iterations, the switching threshold is surpassed and we invoke Newton’s method. Figure 13a shows that the Newton iterations oscillate and cannot converge after 15 iterations (maximum permissible Newton iterations). To demonstrate the difficulty of this problem we revert to SS and continue iterating. An additional 4732 SS iterations are required to converge to a tolerance of < 10−6 , bringing the total SS iteration count to 5742. Figure 13b shows the results using our SS-Newton-TR algorithm. After the maximum number of SS iterations (20 SS steps), the residual is 0.0495. At this point, we switch from SS to Newton’s method. However, the Newton method is infeasible in this situation and we switch to the TR solver immediately. We attain convergence after just 13 TR iterations. We set the default switch criterion to switch = 0.01 for flash calculations, and enforce a maximum upper limit of 20 SS iterations. Continued SS iterations provide little assistance to the TR-solver. For example, in Figure 13a, if we switch SS to TR at switch = 10−5 , we would still require an additional 8 TR iterations for convergence; a modest improvement in the number of TR iterations facilitated by almost 1000 extra SS steps. For a three-phase flash calculation in the critical region, we encounter the convergence problem for Newton’s method because the Hessian matrix is near singular, which is similar to the convergent difficulties in the critical region for a two-phase vapor-liquid system. For example, consider the point (0.893, 85.72) in the JEMA fluid. This point is close to

29



a bi-critical point, and Newton’s method fails for the three-phase flash calculation, unless the residual is first reduced to switch = 2.34 × 10−5 . Figure 14a illustrates these convergent challenges, showing that 1659 SS iterations are required before we can use Newton’s method. The Newton iterations relatively quickly (5 iterations). Using SS alone requires 3649 iterations for convergence. Figure 14b shows a vast improvement in performance, using the SSI-Newton-TR algorithm with the default switch criterion set to switch = 0.01. After just six SS iterations, we switch to the Newton method, which is infeasible. The TR method is then used and we realize convergence in just eight TR iterations. Note that the TR solver advances the solution via SS update at the 8th and 10th iteration in Figure 14b. The SSTR algorithm solves challenging three-phase flash calculations in the critical region without any difficulty. Li and Firoozabadi 29 used SS-Newton for two-phase and three-phase flash calculations. They suggested reverting to SS iterations until convergence should Newton’s method encounter difficulty. We tested this strategy for the same example: xCO2 = 0.893, P = 85.72 bar, for the JEMA fluid. Convergence to the solution requires 3649 SS iterations using the Li and Firoozabadi 29 method. We contend that standalone use of SS in two-phase and three-phase flash calculations in the critical region is prohibitively slow for reservoir simulation. In all our computations, we use a consistent residual threshold for switching from SS to Newton (see Figure 1). For phase-stability testing we use a default value of switch = 0.1, whereas for two-phase and three-phase flash calculations we use a default criterion of switch = 0.01. The selection of these default switching criteria is not arbitrary. We also investigated the effect of switch on execution time. For phase-stability testing, we performed the analysis in the two-phase region for each of the nine fluids (see Figure 6). We varied the switch criterion for the stability tests through 0.01, 0.1 and 1.0. The criterion for the two-phase flash calculations was held constant at switch = 0.01. Taking the average execution time across the nine fluids, we benchmarked performance against the default switch criterion of switch = 0.1. Table 4 lists relative performance. We performed the

30


Page 30 of 56

Page 31 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


same benchmarking for two-phase and three-phase flash calculations. The switch criterion was varied through values of 0.1, 0.01 and 0.001, while the default switch criterion for stability testing was held constant at switch = 0.1. Performance of the two-phase flash was assessed by measuring computation time for data points in the two-phase region. We used the data points in the three-phase region to assess performance of the three-phase flash (see Figure 5). We performed phase equilibrium calculations using several million points across the nine-fluids for each switch criterion. Table 4 lists the results of the analysis. Table 4 shows that phase-stability testing with the default switch criteria of switch = 0.1 results in a 10% decrease in CPU time compared with alternate criteria of switch = 1.0 or switch = 0.01. Total CPU time shows slight sensitivity to the switch value used in two-phase and three-phase flash calculations. In our algorithms, we limit the maximum SS iterations to 20. If no restriction is placed on the maximum number of SS iterations, a small switch value would require a huge number of SS iterations to reach the threshold switch in difficult conditions. Clearly, a strict switch criterion would increase the execution time. We circumvent this difficulty by placing an upper bound on the number of SS iterations. In contrast, if a large switch value were adopted the early switch from SS to Newton would cause some problems for the Newton solver. The implication would be detection of an infeasible condition for Newtons method, and the algorithm would quickly switch to TR iterations. Therefore, the switch value has little effect on the computational performance of both phase stability testing and multiphase flash calculations. The caveat here is that our analysis applies to the SS-Newton-TR (Figure 1). Our conclusion no longer holds for the SS-Newton algorithm without TR. In this circumstance, a much smaller switch value is required to ensure the convergence of the Newton iterations. Table 5 lists the execution times for calculation of the three-phase diagrams shown in Figs. 5a-5i. The calculations were performed using a desktop PC with Intel i7-4790 CPU @3.60GHz processor with 16 GB RAM. The operating system is 64-bit Windows 8.1. The

31



default switch criteria are used for the computations, i.e. switch = 0.1 for stability testing and switch = 0.01 for two-phase and three-phase flash calculations. Table 6 presents information on the total number of three-phase test points for the fluids displayed in Figs. 5a-5i. Table 6 also provides metrics on the iterations used in the threephase region for each test fluid. The default switch criterion switch = 0.01 is used for the three-phase flash calculations. The average number of iterations used for a three-phase flash is 9.2 for SS, 2.4 for Newton and 0.6 for TR, respectively. The TR iterations are higher in ResA and ResB, compared to the other seven fluids. The average number of TR iterations is 0.086 for the seven fluids if we exclude ResA and ResB from Table 6. Table 7 details the results of three-phase flash calculations with the switch criterion set to switch = 0.001 for ReaA and ResB. In these calculations, we remove the limit on the maximum number of SS iterations (previously set to 20). Compared to the results in Table 6, Table 7 shows that the TR iterations decrease significantly while the SS iterations increase enormously. The execution time is 158 seconds for ResA and 165 seconds for ResB. The values are very close to those shown in Table 5. We investigated the effect of the initial trust-region radius ∆0 value on performance and found that the initial value has little effect. The initial value proposed by Petitfrere and Nichita 19 and a simple constant value, for example, ∆0 = 0.25 yield comparable performance in both stability testing and multiphase flash calculations. Values between 0.2 to 0.5 are acceptable. In all our calculations including the solution of the RR equations, we set ∆0 = 0.25.

Solution of Three-phase Rachford-Rice Equations In this section we discuss the efficacy of our Rachford-Rice TR solver. Our TR implementation has proven robust and insensitive to the initial guesses for phase amounts. Table 8 lists the KiV and KiL2 values from the JEMA fluid at point (0.893 85.72) (the same point in the previous three-phase flash test, detailed in Figs. 14a-14b). The KiV values are from 32


Page 32 of 56

Page 33 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


two-phase stability testing and the KiL2 values from two-phase flash. Figure 15a shows the iteration number vs. the residual norm (Euclidean norm of Eq. 15 and Eq. 16) for different initial V and L2 values. In this example, seven iterations are required for convergence if we use a default initial guess of 1/N p = 0.33 for phase amounts. If the initial values are V = 0.33, L2 = −0.33, ten iterations are required. If we set the initial values using the Li and Firoozabadi 29 approach, L2 = 0.83 (from the two-phase flash) and V = 0.01, only one iteration is required. Figs. 15b-15c show that the phase amounts can be negative while the Rachford-Rice equations remain converged. Table 9 lists the KiV values (from two-phase stability testing) and KiL2 values (from two-phase flash) for the MSO fluid at point (0.9483 68.54). The solution of the Rachford-Rice equations is in the two-phase region for this point because the mole fraction of the vapor phase is negative at convergence (see Fig 16b). Fig16a shows the convergence of the TR RR solver for different initial estimates of phase amounts. Again, fewer iterations are required when initial values (V = 0.91 from two-phase flash, L2 = 0.01) proposed by Li and Firoozabadi 29 are used. Figure 16b shows that the V value becomes negative as the iterations proceed, regardless of the initial estimates for phase amounts. Our analysis shows that the starting guess for the phase amounts has a strong effect on the total number of iterations required by our Rachford-Rice TR solver. The phase amounts proposed by Li and Firoozabadi 29 constitute an excellent initial guess for the solution of the three-phase Rachford-Rice problem (Eqs. 15 and 16). However, this approach relies upon the results of two-phase flash and two-phase stability testing calculations. In other words, sequential stability testing and flash calculations starting from a single-phase mixture are required for the first Rachford-Rice iteration (see Figure 4). In reservoir simulation, sequential application of phase equilibrium calculations is always used in the initialization of fluid status in each cell. On the other hand, as the simulation advances in time K-values and phase amounts from the previous time step are used as initial estimates for flash calculations. In this context, the initial estimates for phase amounts proposed by Li and Firoozabadi 29

33



are unavailable. As we show in Figs. 15 and 16, with our TR solver the V and L2 values can be out of physical range (< 0, or > 1) during the iterations. This is important because we frequently encounter negative phase amounts during the iterations when two-phase stability testing provides a false indication of instability. The previous example of the MSO fluid at (0.9483 68.54) (see Table 3) serves as a case in point. An alternative approach to the three-phase RR problem is the use of a bisection method. But a negative 2D bisection method is too slow to be used for the compositional simulation of the three-phase hydrocarbon-CO2 systems. The initial estimate of phase amounts for the solution of the RR equations only affects the number of RR iterations on the first SS iteration in a three-phase flash calculation. In subsequent SS iterations, we estimate phase amounts from the previous SS iteration and the RR equations can converge within three iterations. The starting values for phase fractions in the RR solver are thus identical, other than for the first SS iteration. It follows that starting guesses for phase amounts in the first RR-iteration have no material impact on the cost of three-phase flash calculations.

Conclusions We presented a comprehensive framework on multiphase equilibrium calculations for compositional simulation of CO2 EOR in low temperature reservoirs. Our framework encompasses strategies for reliable stability testing, in addition to prescriptive algorithms for the solution of the TPD and Gibbs free energy minimization problems. We also present a new TR-based solver for robust solution of the Rachford-Rice equations. Herein we provide a brief summary of the key contributions of this paper. We analyzed the shortcomings of current practices in phase stability testing, including: (i) Redundant initial K-value estimates; (ii) incorrectly testing only one phase in two-phase stability testing; (iii) inadequate initial K-value estimates for two-phase stability testing.

34


Page 34 of 56

Page 35 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Prior to this research, consistent detection of three-phase equilibrium for hydrocarbon-CO2 mixtures was not possible. We proposed a new initial estimate for two-phase stability testing, supplementing existing estimates and enabling reliable detection of three-phase behavior. We also significantly reduced the number of initial estimates required; five estimates are adequate for single-phase stability testing and seven are sufficient for two-phase stability testing. We showed that testing of both phases is necessary in two-phase stability testing. Through extensive testing, we demonstrated consistent resolution of two- and three-phase regions throughout pressure-composition space for nine characterized fluids. Complex phase behavior of three-phase hydrocarbon-CO2 systems requires a globally convergent Newton method for both phase stability testing and multiphase flash calculations. We implemented a TR optimization method for the solution of particularly difficult phase equilibrium problems. We demonstrate excellent performance of our implementation for the most challenging problems, including close to the stability test limit locus in stability testing and in the critical region for two-phase and three-phase flash. Our sequential SS-NewtonTR algorithm performs efficiently and robustly in the solution of multiphase equilibrium computations for three-phase hydrocarbon-CO2 systems. Solution of three-phase RR equations is an important element of the three-phase flash calculation. We implemented a trust-region algorithm to solve the RR equations. We demonstrated excellent performance of our RR solver. Performance is reliable and convergence is insensitive to the initial estimates of phase amounts. In addition, our TR solver permits negative values of phase amounts throughout the RR iterations by minimizing a convex function with constraints. We have tested tens of millions of points across single-phase, two-phase, and three-phase regions with a test suite comprising nine fluids from the literature. We have not encountered a single failure throughout our rigorous testing. Our results demonstrate that the framework in this research is applicable to multiphase equilibrium calculations in the compositional simulation of CO2 injection in low temperature reservoirs.

35



Acknowledgement We acknowledge the financial support from the SUPRI-B industrial affiliates program at Stanford University.

References (1) Brock, W.; Bryan, L. Summary results of CO2 EOR field tests, 1972-1987. low permeability reservoirs symposium. 1989. (2) Godec, M. L.; Kuuskraa, V. A.; Dipietro, P. Opportunities for using anthropogenic CO2 for enhanced oil recovery and CO2 storage. Energy & Fuels 2013, 27, 4183–4189. (3) Melzer, L. S. Carbon dioxide enhanced oil recovery (CO2 EOR): Factors involved in adding carbon capture, utilization and storage (CCUS) to enhanced oil recovery. Center for Climate and Energy Solutions 2012, (4) Nghiem, L.; Li, Y. Effect of phase behavior on CO2 displacement efficiency at low temperatures: model studies with an equation of state. SPE Reservoir Engineering 1986, 1, 414–422. (5) Pan, H.; Chen, Y.; Sheffield, J.; Chang, Y.-B.; Zhou, D. Phase-Behavior Modeling and Flow Simulation for Low-Temperature CO2 Injection. SPE Reservoir Evaluation & Engineering 2015, 18, 250–263. (6) Okuno, R.; Johns, R. T.; Sepehrnoori, K. Three-phase flash in compositional simulation using a reduced method. SPE Journal 2010, 15, 689–703. (7) Michelsen, M. L. The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilibria 1982, 9, 21–40. (8) Michelsen, M. Phase equilibrium calculations. What is easy and what is difficult? Computers & chemical engineering 1993, 17, 431–439. 36


Page 36 of 56

Page 37 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(9) Hoteit, H.; Firoozabadi, A. Simple phase stability-testing algorithm in the reduction method. AIChE journal 2006, 52, 2909–2920. (10) Nichita, D. V.; Broseta, D.; Montel, F. Calculation of convergence pressure/temperature and stability test limit loci of mixtures with cubic equations of state. Fluid Phase Equilibria 2007, 261, 176–184. (11) Murray, W. Second derivative methods. Numerical Methods for Unconstrained Optimization (Ed. W. Murray). Academic Press, London 1972, (12) Trangenstein, J. A. Customized minimization techniques for phase equilibrium computations in reservoir simulation. Chemical Engineering Science 1987, 42, 2847–2863. (13) Perschke, D. R. Equation-of-state phase-behavior modeling for compositional simulation. Ph.D. thesis, Texas Univ., Austin, TX (USA), 1988. (14) Jindrová, T.; Mikyˇska, J. General algorithm for multiphase equilibria calculation at given volume, temperature, and moles. Fluid Phase Equilibria 2015, 393, 7–25. (15) Smejkal, T.; Mikyˇska, J. Phase stability testing and phase equilibrium calculation at specified internal energy, volume, and moles. Fluid Phase Equilibria 2017, 431, 82–96. (16) Nichita, D. V. Fast and robust phase stability testing at isothermal-isochoric conditions. Fluid Phase Equilibria 2017, 447, 107–124. (17) Nichita, D. V. Volume-based phase stability testing at pressure and temperature specifications. Fluid Phase Equilibria 2018, 458, 123–141. (18) Nichita, D. V. A volume-based approach to phase equilibrium calculations at pressure and temperature specifications. Fluid Phase Equilibria 2018, 461, 70–83. (19) Petitfrere, M.; Nichita, D. V. Robust and efficient trust-region based stability analysis and multiphase flash calculations. Fluid Phase Equilibria 2014, 362, 51–68. 37



(20) Deuflhard, P. Newton methods for nonlinear problems: affine invariance and adaptive algorithms; Springer Science & Business Media, 2011; Vol. 35. (21) Nghiem, L. X.; Aziz, K.; Li, Y. A robust iterative method for flash calculations using the Soave-Redlich-Kwong or the Peng-Robinson equation of state. Society of Petroleum Engineers Journal 1983, 23, 521–530. (22) Lucia, A.; Liu, D. An acceleration method for dogleg methods in simple singular regions. Industrial & engineering chemistry research 1998, 37, 1358–1363. (23) Powell, M. J. An efficient method for finding the minimum of a function of several variables without calculating derivatives. The computer journal 1964, 7, 155–162. (24) Moré, J. J.; Sorensen, D. C. Computing a trust region step. SIAM Journal on Scientific and Statistical Computing 1983, 4, 553–572. (25) Pan, H.; Firoozabadi, A. Complex multiphase equilibrium calculations by direct minimization of Gibbs free energy by use of simulated annealing. SPE Reservoir Evaluation & Engineering 1998, 1, 36–42. (26) McKinnon, K.; Mongeau, M. A generic global optimization algorithm for the chemical and phase equilibrium problem. Journal of Global Optimization 1998, 12, 325–351. (27) Nghiem, L. X.; Li, Y.-K. Computation of multiphase equilibrium phenomena with an equation of state. Fluid Phase Equilibria 1984, 17, 77–95. (28) Computer Modelling Group Ltd, WinProp User Guide. 2011. (29) Li, Z.; Firoozabadi, A. General strategy for stability testing and phase-split calculation in two and three phases. SPE Journal 2012, 17, 1096–1107. (30) Nichita, D. V.; Petitfrere, M. Phase equilibrium calculations with quasi-Newton methods. Fluid Phase Equilibria 2015, 406, 194–208. 38


Page 38 of 56

Page 39 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(31) Haugen, K. B.; Firoozabadi, A.; Sun, L. Efficient and robust three-phase split computations. AIChE Journal 2011, 57, 2555–2565. (32) Firoozabadi, A.; Pan, H. Fast and Robust Algorithm for Compositional Modeling: Part I-Stability Analysis Testing. SPE Journal 2002, 7, 78–89. (33) Pan, H.; Firoozabadi, A. Fast and Robust Algorithm for Compositional Modeling: Part II-Two-Phase Flash Computations. SPE Journal 2003, 8, 380–391. (34) Li, Y.; Johns, R. T. Rapid flash calculations for compositional simulation. SPE Reservoir Evaluation & Engineering 2006, 9, 521–529. (35) Petitfrere, M.; Nichita, D. V. Multiphase equilibrium calculations using a reduction method. Fluid Phase Equilibria 2015, 401, 110–126. (36) Rachford Jr, H.; Rice, J. Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium. Pet. Trans. AIME 1952, 195, 237–238. (37) Okuno, R. Modeling of multiphase behavior for gas flooding simulation. Ph.D. thesis, The University of Texas at Austin, 2009. (38) Leibovici, C. F.; Neoschil, J. A solution of Rachford-Rice equations for multiphase systems. Fluid Phase Equilibria 1995, 112, 217–221. (39) Michelsen, M. Calculation of multiphase equilibrium. Computers & chemical engineering 1994, 18, 545–550. (40) Leibovici, C. F.; Nichita, D. V. A new look at multiphase Rachford–Rice equations for negative flashes. Fluid Phase Equilibria 2008, 267, 127–132. (41) Yan, W.; Stenby, E. H. On multiphase negative flash for ideal solutions. Fluid Phase Equilibria 2012, 322, 41–47.

39



(42) Li, Z.; Firoozabadi, A. Initialization of phase fractions in Rachford–Rice equations for robust and efficient three-phase split calculation. Fluid Phase Equilibria 2012, 332, 21–27. (43) Iranshahr, A.; Voskov, D.; Tchelepi, H. Generalized negative-flash method for multiphase multicomponent systems. Fluid Phase Equilibria 2010, 299, 272–284. (44) Okuno, R.; Johns, R.; Sepehrnoori, K. A new algorithm for Rachford-Rice for multiphase compositional simulation. SPE Journal 2010, 15, 313–325. (45) Wilson, G. M. A modified Redlich-Kwong equation of state, application to general physical data calculations. 65th National AIChE Meeting, Cleveland, OH. 1969. (46) Michelsen, M. L. The isothermal flash problem. Part I. Stability. Fluid phase equilibria 1982, 9, 1–19. (47) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs energy analysis of phase equilibria. Society of Petroleum Engineers Journal 1982, 22, 731–742. (48) Zhou, Y. Parallel General-Purpose Reservoir Simulation with Coupled Reservoirs and Multisegment Wells. Ph.D. thesis, Stanford Univ., Stanford, CA (USA), 2012. (49) Robinson, D. B.; Peng, D.-Y. The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs; Gas Processors Association, 1978. (50) Michelsen, M.; Mollerup, J. Thermodynamic Models: Fundamental & Computational Aspects; Tie-Line Publications, 2004. (51) Nocedal, J.; Wright, S. J. Sequential quadratic programming; Springer, 2006. (52) Pan, H.; Imai, M.; Connolly, M.; Tchelepi, H. Efficient and robust solver for multiphase Rachford-Rice equations in compositional and thermal simulations. SPE Reservoir Simulation Conference Symposium. 2019. 40


Page 40 of 56

Page 41 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(53) Connolly, M. An isenthalpic-based compositional framework for nonlinear thermal simulation. Ph.D. thesis, Stanford Univ., Stanford, CA (USA), 2018. (54) Khan, S.; Pope, G.; Sepehrnoori, K. Fluid characterization of three-phase CO2 /oil mixtures. SPE/DOE Enhanced Oil Recovery Symposium. 1992.

41



Page 42 of 56

Tables Table 1: Three-phase Points only Detected by Initial Estimate KiCO2 −95 JEMA xCO2 P (bar) 0.703 83.82 0.699 82.16 0.727 81.56 0.754 81.44 0.766 81.42

MRO xCO2 P (bar) 0.658 90.84 0.656 91.24 0.646 93.4 0.654 92.08 0.654 91.92 0.654 91.72

NWE xCO2 P (bar) 0.572 91.16 0.568 91.68 0.566 92.12 0.566 92.16

ResB xCO2 P (bar) 0.675 89.87 0.836 86.96 0.837 86.94 0.838 86.92 0.84 86.87 0.841 86.84

Table 2: Three-phase Points only Detected by the Second Two-phase Flash Result NWE xCO2 P (bar) 0.793 69.46 0.795 69.34 0.8 69.1

MSO xCO2 P (bar) 0.9483 68.54 0.95 68.57 0.9553 68.7

OilB xCO2 P (bar) 0.854 75.08 0.855 75

Table 3: Calculation Results at xCO2 = 0.793, P = 69.46 bar for NWE TPD (1P stability) -0.28819 -0.07871

G(2P flash) -2.81024 -2.81542

TPD (2P stability) -0.07871 -0.0009

Table 4: Effect of switch on Execution Time switch 1.0 0.1 0.01 0.001

Stability 2P-flash 1.10 1.00 1.04 1.08 1.00 1.02

42

3P-flash 1.00 1.00 0.983


G(3P flash) failed -2.81544

Page 43 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5: Execution Time (seconds) in the Computation of the Three-phase P-x Diagrams Shown in Figs. 5a-5i Fluid Ntest Time

BSB 338701 255

JEMA 203401 163

NWE 630450 413

OilG 625500 437

MRO 500400 740

MSO 360400 208

OilB 255300 727

ResA 300300 155

ResB 353000 160

Ntest: total number of calculated points. Table 6: Total and Average Iterations in the Computation of the Three-phase Region Shown in Figs. 5a-5i Fluid

N3P

BSB JEMA MSO MRO NEW OilB OilG ResA ResB average

54636 50244 31658 78215 97806 29570 104241 23330 22173

Total Iterations SSI Newton TR 537205 144814 12148 555085 128743 6726 237663 60881 172 887844 181704 9888 867252 220142 5603 234259 64846 917 1078506 232387 2524 190167 65544 40392 173230 57561 68239

Average Iterations SSI Newton TR 9.83 2.65 0.222 11.05 2.56 0.134 7.51 1.92 0.005 11.35 2.32 0.126 8.87 2.25 0.057 7.92 2.19 0.031 10.35 2.23 0.024 8.15 2.81 1.731 7.81 2.6 3.078 9.2 2.4 0.6

N3P: total number of calculated three-phase flash. The default switch criterion switch = 0.01 from SS to Newton is used. Table 7: Total and Average Iterations in the Computation of the Three-phase Region Shown in Figs. 5h and 5i Fluid

N3P

ResA ResB

23330 22173

Total iterations SSI Newton TR 778529 54935 2861 1172888 54722 3614

Average Iterations SSI Newton TR 33.37 2.35 0.123 52.90 2.47 0.163

N3P: total number of calculated three-phase flash. The switch criterion switch = 0.001 from SS to Newton is used. Table 8: Phase Equilibrium Constants for RR Equations from JEMA KiL2 1.544580305 2.336017325 8.869147335E-01 3.712005301E-01 4.768114520E-02 1.059266958E-03 1.766173590E-06

KiV 1.546336826 2.379576919 8.818093562E-01 3.597420492E-01 4.404202880E-02 8.924834994E-04 1.302131664E-06 43



Table 9: Phase Equilibrium Constants for RR Equations from MSO KiL2 1.479781968 6.117039468E-02 1.625155948E-02 2.749741486E-03 2.664184498E-04 1.309992275E-05 1.928151262E-08

KiV 1.420340741 2.964408254E-01 1.805854981E-01 9.303281846E-02 3.891826286E-02 1.263057652E-02 1.106068886E-03

44


Page 44 of 56

Page 45 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figures

Figure 1: Schematic illustration of the sequential application of SS, Newton and trust region solvers. Numerical solution of both phase stability testing and flash calculations proceeds in this procedure. The threshold for convergence of the residual norm is denoted conv . The term switch refers to the threshold residual at which the SS iterations cease, and the algorithm switches to the Newton method .

Figure 2: Algorithmic description for updating the trust-region radius ∆k .

45



Figure 3: Algorithmic description for the solution of the TR sub-problem. λ1 is the smallest eigenvalue of the Hessian matrix H.

Figure 4: Schematic illustration of the solution strategy for multiphase equilibrium calculations. First, we perform two-phase flash with the K-values that correspond to the lowest negative TPD in single-phase stability testing. This is depicted on the left. If the initial three-phase flash fails, another two-phase flash is performed. We initiate the second twophase flash using the K-values associated with the next lowest negative TPD in single-phase stability testing. Two-phase stability testing and three-phase flash (if required) then follows.

46


Page 46 of 56

Page 47 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a) BSB

(b) JEMA

(c) NWE

(d) OilG

(e) MRO

(f) MSO

(g) OilB

(h) ResA

(i) ResB

Figure 5: P-x diagrams focused on the three-phase region for nine test fluids. The terms 1P, 2P and 3P represent single-phase, two-phase and three-phase regions, respectively. The resolution is ∆P = 0.01 or 0.02 bar, and ∆x = 0.001 in these figures.

47



(a) NWE

(b) OilG

Page 48 of 56

(c) ResA

Figure 6: P-x diagrams focused on the two-phase liquid-liquid (L-L2 ) region for (a) NWE, (b) OilG and (c) ResA. The terms 1P, 2P and 3P are used to represent the single-phase, twophase and three-phase regions, respectively. The resolution of the P-x diagrams is ∆P = 1 bar, and ∆x = 0.001. A single-phase region is present at moderate concentrations of injected CO2 . A very narrow single-phase region is present on the far right of each phase diagram at high concentrations of injected CO2 . The single-phase that appears in this region is a CO2 -rich liquid phase (L2 ).

(a) NWE

(b) OilG

(c) MSO

Figure 7: P-x diagram focused on the region of high CO2 concentration for (a) NWE, (b) OilG and (c) MSO. The terms 1P, 2P and 3P represent single-phase, two-phase and three-phase regions, respectively. The resolution is ∆P = 1 bar, ∆x = 0.0001 in these P-x diagrams. The CO2 -rich liquid phase (L2 ) is present on the far right of each diagram, corresponding to very high concentrations of CO2 .

48


Page 49 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a) Test L Phase

(b) Test V Phase

(c) Test V and L Phases

Figure 8: P-x diagrams constructed using Nc+4 initial estimates for two-phase stability testing of the BSB fluid, when the test phase is taken to be (a) the heavy L-phase; (b) the light V-phase; (c) both the L and V-phase. From these P-x diagrams it is clear that two-phase stability testing fails to identify points in the 3-phase region when only one of the two phases is used as a test-phase. In case (a) when only the L-phase is tested, three-phase behavior is not identified at moderate concentrations of CO2 . When the V-phase serves as the test phase in case (b), some points in 3-phase region are incorrectly identified as twophase at high CO2 concentrations close to bi-critical point. In addition, a few points are missing between x = 0.65 − 0.8. In case (c) the entire three-phase region is resolved when both the V and L-phases are used as test phases.

(a) Test L Phase

(b) Test V Phase


Figure 9: P-x diagrams for JEMA fluid, illustrating the importance of testing both phases in two-phase stability testing. The phase diagrams were constructed with Nc+4 initial estimates used for two-phase stability testing. Inconsistent identification of three-phase behavior is clearly visible for the P-x phase diagram. In case (a) simply testing the heavy L-phase is insufficient. In case (b) only the light V-phase is tested for instability. The testing fails to detect 3-behavior occurs at high concentrations of injected CO2 , close to bi-critical point, and at a few points between x = 0.7 and x = 0.8. For case (c) almost all of the three-phase region is resolved when both of the phases are tested for stability. But a few points are still incorrectly identified as 2-phase between x = 0.65 and x = 0.8 in the 3-phase region.

49



(a) Test L Phase

(b) Test V Phase


Figure 10: P-x diagrams for NWE fluid. The diagrams were constructed using Nc+4 initial estimates in two-phase stability testing. In case (a) the L-phase serves as the test phase for the two-phase stability tests. Many points in the 3-phase region are incorrectly identified as two-phase for moderate concentrations of injected CO2 . In case (b) the V-phase is the test-phase. Points in the 3-phase region are incorrectly identified as 2-phase, this time at x < 0.6 and between x = 0.7 and 0.85. Finally, in case (c) two-phase stability testing using both the L and V-phase misses some points in the 3-phase region at x < 0.65.

(b) ResB

(a) MRO

Figure 11: P-x diagrams for (a) MRO, and (b) ResB. These P-x diagrams were constructed using a total of Nc+4 initial estimates for two-phase stability testing. Here, we only show the results of two-phase stability testing using both the L and V phases as the test phase. In (a) many points in the 3-phase region are incorrectly identified as 2-phase when the injected CO2 concentration is between 0.6 and 0.7. For the ResB fluid in case (b), 3-phase equilibrium is not detected by the two-phase stability testing between x = 0.8 and 0.9.

50


Page 50 of 56

Page 51 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(b) SS-Newton-TR

(a) SS-Newton

Figure 12: Number of iterations for stability testing calculation with XCO2 = 0.603, P = 82.4 bar for OilG. In case (a) the SS-Newton method is used, SS iterations reduce the magnitude of the residual norm very slowly and the Newton iterations diverge when the switch criterion is switch = 10−5 . In case (b) the SS-Newton-TR algorithm is used. After two SS iterations, the TR solver converges in 11 iterations.

(a) SS-Newton

(b) SS-Newton-TR

Figure 13: Number of iterations for two-phase flash calculation at XCO2 = 0.747, P = 77.06 bar for OilG. In case (a) the SS-Newton method is used, the SS iterations reduce the residual norm very slowly, and the Newton iterations diverge, even with a stringent switching criterion of switch = 10−5 . In case (b) the SS-Newton-TR algorithm is used. After 20 SS iterations, TR converges in just 13 iterations.

51



(b) SS-Newton-TR

(a) SS-Newton

Figure 14: Number of iterations for three-phase flash calculation at XCO2 = 0.893, P = 85.72 bar for JEMA. In case (a) the SS-Newton method is used, the SS iterations reduce residual slowly and Newton iterations converge only if the switching criterion is switch ≤ 2.34 × 10−5 . In case (b) the SS-Newton-TR algorithm is used. After 6 SS iterations, the TR solver converges in just 8 iterations.

(a) Number of RR iterations for convergence

(b) Vapor phase amount V

(c) Light liquid phase amount L2

Figure 15: Results for solving RR equations using different initial estimates for JEMA. The K-values are listed in Table 8. The number of iterations are sensitive to the initial values. 52


Page 52 of 56

Page 53 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a) Number of RR iterations for convergence

(b) Vapor phase amount V

Figure 16: Results for solving RR equations using different initial estimates for MSO. The K-values are listed in Table 9.

53



Page 54 of 56

Appendix A: Thermodynamic Properties of Fluids Table A1: Component Data of BSB Oil Component CO2 C1 C2−3 C4−6 C7−15 C16−27 C28+

z 0.0337 0.0861 0.1503 0.1671 0.3304 0.1611 0.0713

Tc(K) 304.2 160 344.2 463.22 605.75 751.02 942.48

Pc(bar) 73.76 46 44.99 34 21.75 16.54 16.42

ω 0.225 0.008 0.1305 0.2407 0.6177 0.9566 1.2683

MW(g/mol) 44.01 16.043 37.2002 69.4984 140.956 280.9914 519.6219

ki,CO2 0.055 0.055 0.055 0.105 0.105 0.105

Table A2: Component Data of JEMA Oil Component CO2 C1 C2−3 C4−6 C7−16 C17−29 C30+

z 0.0192 0.0693 0.1742 0.1944 0.3138 0.1549 0.0742

Tc(K) 304.2 166.67 338.81 466.12 611.11 777.78 972.22

Pc(bar) 73.76447 46.00147 45.5271 33.67863 20.95029 15.88218 15.84339

ω 0.225 0.008 0.126 0.2439 0.6386 1.0002 1.2812

MW(g/mol) 44.01 16.043 36.0126 70.5203 147.182 301.4762 562.8054

ki,CO2 0.05 0.05 0.05 0.09 0.09 0.09

Table A3: Component Data of NWE (North Ward Estes) Oil Component CO2 C1 C2−3 C4−6 C7−14 C15−24 C25+

z 0.0077 0.2025 0.118 0.1484 0.2863 0.149 0.0881

Tc(K) 304.2 190.6 343.64 466.41 603.07 733.79 923.2

Pc(bar) ω 73.76 0.225 46 0.008 45.05 0.1301 33.5 0.2436 24.24 0.6 18.03 0.903 17.26 1.229

54


MW(g/mol) 44.01 16.043 38.3985 72.824 135.8191 257.7499 479.9548

ki,CO2 0.12 0.12 0.12 0.09 0.09 0.09

Page 55 of 56 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table A4: Component Data of OilG Component CO2 C1 C2−3 C4−6 C7−14 C15−25 C26+

z 0.0169 0.1752 0.2244 0.1673 0.2422 0.1216 0.0524

Tc(K) 304.2 174.44 347.26 459.74 595.14 729.98 910.18

Pc(bar) 73.76 46 44.69 34.18 21.87 16.04 15.21

ω 0.225 0.008 0.1331 0.2358 0.5977 0.9118 1.2444

MW(g/mol) 44.01 16.043 37.9086 68.6715 135.0933 261.103 479.6983

ki,CO2 0.085 0.085 0.085 0.104 0.104 0.104

Table A5: Component Data of MRO (Maljamar Reservoir Oil) Component CO2 C1 C2 C3 nC4 C5−7 C8−10 C11−14 C15−20 C21−28 C29+

z

Tc(K) 304.211 190.6 305.4 369.8 425.2 516.667 590 668.611 745.778 812.667 914.889

0 0.2939 0.1019 0.0835 0.0331 0.1204 0.1581 0.0823 0.0528 0.0276 0.0464

Pc(bar) 73.819 45.4 48.2 41.9 37.5 28.82 23.743 18.589 14.8 11.954 8.523

ω 0.225 0.008 0.098 0.152 0.193 0.2651 0.3644 0.4987 0.6606 0.8771 1.2789

MW(g/mol) 44 16 30.1 44.1 58.1 89.9 125.7 174.4 240.3 336.1 536.7

ki,CO2

ki,C1

0.115 0.115 0.115 0.115 0.115 0.115 0.115 0.115 0.115 0.115

0.045 0.055 0.055 0.06 0.08 0.28

Table A6: Component Data of MSO (Maljamar Separator Oil) Component CO2 C5−7 C8−10 C11−14 C15−20 C21−28 C29+

z 0 0.2354 0.3295 0.1713 0.1099 0.0574 0.0965

Tc(K) 304.211 516.667 590 668.611 745.778 812.667 914.889

Pc(bar) 73.819 28.82 23.743 18.589 14.8 11.954 8.523

55

ω 0.225 0.2651 0.3644 0.4987 0.6606 0.8771 1.2789


MW(g/mol) 44 89.9 125.7 174.4 240.3 336.1 536.7

ki,CO2 0.115 0.115 0.115 0.115 0.115 0.115


Page 56 of 56

Table A7: Component Data of OilB Component CO2 N2 C1 C2 C3 IC4 N C4 IC5 N C5 C6 P C1 P C2 P C3 P C4 P C5 P C6

z 0.0011 0.0048 0.163 0.0403 0.0297 0.0036 0.0329 0.0158 0.0215 0.0332 0.184152 0.163891 0.127256 0.096888 0.058808 0.023105

Tc(K) 304.2 126.2 190.6 305.4 369.8 408.1 425.2 460.4 469.6 506.35 566.55 647.06 719.44 784.93 846.33 919.39

Pc(bar) 72.8 33.5 45.4 48.2 41.9 36 37.5 33.4 33.3 33.9 25.3 19.1 14.2 10.5 7.5 4.76

ω 0.225 0.04 0.008 0.098 0.152 0.176 0.193 0.227 0.251 0.299 0.3884 0.5289 0.6911 0.8782 1.1009 1.4478

MW(g/mol) 44.01 28.01 16.04 30.07 44.1 58.12 58.12 72.15 72.15 84 112.8 161.2 223.2 304.4 417.5 636.8

ki,CO2

ki,N2

ki,C1

-0.02 0.075 0.08 0.08 0.085 0.085 0.085 0.085 0.095 0.095 0.095 0.095 0.095 0.095 0.095

0.08 0.07 0.07 0.06 0.06 0.06 0.06 0.05 0.1 0.12 0.12 0.12 0.12 0.12

0.003 0.01 0.018 0.018 0.025 0.026 0.036 0.049 0.073 0.098 0.124 0.149 0.181

Table A8: Component Data of ResA Component CO2 C1 C2−3 C4−6 C11 C31

z 0.0203 0.1419 0.1963 0.1742 0.3764 0.0909

Tc(K) 304.19 190.58 337.571 456.435 648.651 898.714

Pc(bar) 73.815 46.043 45.651 35.753 23.482 11.731

MW(g/mol) 44.01 16.043 37.072 68.062 148.3568 424.0175

ω 0.228 0.011 0.125 0.222 0.445 1.034

ki,CO2

ki,C1

ki,C2−3

0.1 0.13 0.125 0.1 0.1

0.005 0.025 0.078 0.098

0.01 0.01 0.01

Table A9: Component Data of ResA Component CO2 C1 C2−3 C4−6 C12 C35

z 0.014 0.073 0.176 0.194 0.402 0.141

Tc(K) 304.19 190.58 347.37 460.27 669.3 915.42

Pc(bar) 73.815 46.043 55.063 34.998 20.465 9.87

ω 0.228 0.011 0.117 0.227 0.504 1.179

MW(g/mol) 44.01 16.043 36.83 69.91 165.9 484.76

ki,CO2

ki,C1

0.1 0.13 0.06 0.09 0.115

0.067 0.083 0.1462 0.1642

ki,C2−3

0.01 0.01 0.01

Table A10: Reservoir Temperature of Test Fluids Fluid T (K)

BSB 313.71

JEMA 316.48

NWE 301.48

OilG 307.59

MRO 305.35 56

MSO 305.35


OilB 307.6

ResA 311.483

ResB 316.483

Multiphase Equilibrium Calculation Framework for ... - ACS Publications

Recommend Documents