Three-Phase Equilibrium Computations for HydrocarbonâWater

Article pubs.acs.org/IECR

Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Three-Phase Equilibrium Computations for Hydrocarbon−Water Mixtures Using a Reduced Variables Method Michael Connolly,*,†,‡ Huanquan Pan,‡ and Hamdi Tchelepi‡ †

McKinsey & Company, Houston, Texas 77002, United States Department of Energy Resources Engineering, Stanford University, Stanford, California 94305, United States

Downloaded via UNIV OF LOUISIANA AT LAFAYETTE on August 5, 2019 at 20:32:26 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: Thermal compositional simulation requires phase-equilibrium calculations for at least three fluid phases. The use of precomputed equilibrium ratios (K-values) has long been justified on the basis of efficiency. However, this method may not appropriately represent thermal recovery mechanisms. An equation of state (EOS) approach is more rigorous, though prohibitively costly. In thermal recovery processes the injection of steam results in hydrocarbon−water interactions at elevated temperatures. Representation of phase behavior for these systems in a nonisothermal context gives rise to a set of challenges owing to the polarity of the water molecule and resulting nonideal phase behavior. In this research we address two difficulties pertinent to thermal compositional reservoir simulation: (i) robust phase-stability analysis using a predetermined set of initial estimates for phase-equilibrium ratios; (ii) numerical solution of phase-split (flash) calculations in the presence of trace components. Phase-stability testing can be difficult for hydrocarbon−water mixtures. Three-phase regions can be extremely narrow in pressure−temperature space, making the resolution of phase boundaries problematic. The standard approach to phase stability testing for hydrocarbon reservoir fluids entails using a series of initial guesses for the equilibrium ratios based on variations of the Wilson correlation (Wilson, G. M. A modif ied Redlich-Kwong equation of state, application to general physical data calculations; 65th National AIChE Meeting, Cleveland, OH, 1969) and trial-phase compositions dominated by each of the Nc components present. For hydrocarbon−water mixtures, fewer initial guesses are required. We propose a physics-based strategy, sensitive to the distinct behavior of water. We use the steam saturation pressure to guide our selection of trial phase compositions. Below saturated steam pressure, only two sets of equilibrium ratios are required to identify the correct phase state. Above the saturation pressure, we expand our set of initial K-value estimates to account for the appearance of a near pure aqueous phase. The K-values obtained from stability analysis are used to perform two-phase flash computations. Then, the stability of the two-phase mixture is assessed. The phase state of the system following the two-phase flash guides the choice of the next set of K-values. This strategy has two direct benefits. First, the number of trial phases required for stability testing is predetermined. Second, the reliability of the phase-stability tests and the ensuing equilibrium calculations is greatly improved. The solution procedure for phase-split calculations in isothermal compositional simulation is typically a two-stage process that uses successive-substitution in conjunction with the Newton method. However, in a nonisothermal setting the presence of trace amounts of hydrocarbons in the aqueous phase can produce ill-conditioned linearized systems when standard (conventional) variables are used in the Newton loop. We have developed a reduced order model for three-phase split calculations that allows us to circumvent these numerical difficulties. The reduced method is formulated to take advantage of the sparsity of the binary interaction parameter matrix. For the first time, we demonstrate the efficacy of a reduced-variables formulation for phase-split calculations involving trace components. We present our protocol for hydrocarbon−water phase-equilibrium computations through comprehensive testing of characterized fluids from the literature. Through our approach to robust phase-stability testing and phase-split calculations, we are able to consistently resolve phase boundaries.

■

INTRODUCTION

involve complex multiphase, multicomponent interactions in porous media. The phase behavior of these systems involves water and hydrocarbon components that partition across multiple fluid phases as a function of composition, pressure,

Over the last 40 years much of the research into efficient and robust phase-equilibrium calculations has been motivated by the study of miscible recovery mechanisms, where CO2 or other gaseous injectants are used to achieve miscibility and improve displacement efficiency. Increasingly, there is an emphasis on EOS representation of phase equilibria for other mechanisms of enhanced oil recovery, including chemical2 and thermal recovery processes.3,4 Thermal recovery processes © XXXX American Chemical Society

Received: Revised: Accepted: Published: A

February 3, 2019 July 13, 2019 July 22, 2019 July 22, 2019 DOI: 10.1021/acs.iecr.9b00695 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research and temperature. The widely accepted industry practice is to use precomputed equilibrium ratios (K-values) that are independent of composition to represent the underlying phase behavior of hydrocarbon−water systems in thermal settings. Compositional independence is a reasonable assumption in the presence of linear phase boundaries.5 However, recent EOS-based thermal compositional simulators have shown that the K-values representation of phase behavior can result in substantial errors in simulated oil recovery.6−8 For this reason, there is a need for EOS-based thermal compositional simulation. In most miscible EOR processes, miscibility is achieved in a multicontact displacement in which compositional routes traverse the critical locus.9 As a result, phase-equilibrium computations have been developed with a view to resolving phase behavior in the near-critical region.10 Although water is always present in petroleum reservoirs and in contact with hydrocarbon phases, its presence has little impact on the phase behavior of hydrocarbon fluids in an isothermal setting. For this reason, water is typically treated as an inanimate phase in compositional reservoir simulation. In contrast, water is the thermodynamically dominant component in thermal EOR processes such as steam-flooding and steam-assisted gravity drainage (SAGD), driving viscosity reduction of oil through steam latent heat transfer. Phase transitions in thermal recovery occur frequently and abruptly, owing to both the sharp temperature gradients and narrow-boiling-point behavior of hydrocarbon−water mixtures. Rigorous phase-behavior modeling is an important element of thermal compositional simulation when studying complex recovery mechanisms.8 In reservoir simulation, hydrocarbon phase equilibrium is represented using cubic equations of state such as the Peng− Robinson (PR) EOS11,12 and Soave−Redlich−Kwong (SRK) EOS.13 Cubic EOS phase-equilibrium computations are welldeveloped for nonpolar hydrocarbon mixtures. The representation of hydrocarbon−water mixtures using an EOS remains challenging due to the polarity of the water molecule. The strong hydrogen bonding gives rise to distinctive, nonideal behavior. Several authors have focused on modifications to cubic EOS models for improved representation of hydrocarbon−water systems.14,15 In order of increasing complexity, these modifications include the use of dual sets of binary interaction coefficients, asymmetric mixing rules, and cubicplus equations of state. Changes to the underlying cubic EOS provide for more accurate representation of phase behavior and have found some limited application in research compositional simulators.16 However, in practice they are seldom implemented in commercial reservoir simulation models. In this paper we focus our attention on algorithms for hydrocarbon−water phase-equilibrium computations using the PR EOS without modification to the van Der Waals mixing rules. The flash problem involves two constitutive uncertainties: (i) the number of phases present; (ii) the composition and amount of each phase present. The first pertains to stability analysis, and the second embodies the phase-split problem. A key challenge in phase-equilibrium calculations is that the number of phases at thermodynamic equilibrium is unknown a priori. There are two broad approaches to addressing this problem:

II. Perform stability-testing and phase-equilibrium calculations in a stepwise procedure. Approaches I and II are both widely used in reservoir simulation. Preallocating the number of phases can be advantageous if assumptions can be made as to the phase state from information at previous time steps or based on the distinctive behavior of the mixture. If an invalid solution is obtained, the number of phases can be decreased until a physical solution is attained. In the case of hydrocarbon−water mixtures, the ubiquitous aqueous phase generally prompts preallocation of the phase state. Flash calculations are generally initiated from a two-phase liquid hydrocarbon−aqueous state. However, this can be a precarious procedure. It may suffer from poor initial estimates of equilibrium ratios and be particularly problematic when phase changes are abrupt and rapid, as in the case of thermal reservoir simulation. In this paper we show that it may be difficult to identify three-phase behavior following the assumption of a liquid−liquid state. As an alternative, stepwise stability testing provides confidence as to the number of phases, as well as good initial estimates for subsequent phase-split calculations. In addition, stability tests are necessary for initialization of compositional reservoir simulation models when knowledge of the phase state from previous time steps is unavailable. Phase-stability analysis is formulated in terms of a tangent plane distance (TPD) function.17 Stability-testing and phasesplit calculations can be solved via two broad approaches: (i) direct solution of nonlinear equations; (ii) minimization of the TPD function or the Gibbs free energy. Successive substitution is almost always used to initiate the solution process, and the initial iterations on the nonlinear equations usually exhibit rapid convergence. Converging to a strict tolerance is more challenging, and a Newton procedure for option (i) or (ii) is usually required. Global optimization methods have been investigated for phase-stability testing and phase-split calculations.18−21 However, global minimization techniques remain computationally expensive and have not found application in reservoir simulation software. One of the key difficulties with stability-analysis and phasesplit algorithms is the initiation of the calculations. The equilibrium ratios from stability testing can be used as an initial estimate in the phase-split calculation. For stability testing, the widely used industry practice is to perform multiple stability tests using several sets of initial K-value estimates. Practically, convergence is to one of the stationary points of the TPD function. The stationary points, however, may correspond to trivial solutions. To successfully detect instability, heuristic strategies are used to initiate the calculation with a trial phase composition. For hydrocarbon fluids, the Wilson correlation1 has been used extensively to generate initial K-value estimates: ln(K iWilson) = ln

Tc, i yz i zz + 5.373(1 + ωi)jjjj1 − P T z{ k

Pc, i

(1)

where Tc,i and Pc,i refer to the critical temperature and pressure of component i and ω is the acentric factor. Michelsen17 suggested that the trial phase may be assumed to be a pure substance. Li and Firoozabadi22 suggested using multiple estimates derived from the Wilson correlation, in addition to corner point estimates:

I. Preallocate the number of phases and proceed directly to the phase-split calculation.

{K iinitial} = {K iWilson , 1/K iWilson ,

3

K iWilson , 1/ 3 K iWilson , K ic} (2)

B

DOI: 10.1021/acs.iecr.9b00695 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research where Kci refers to Nc initial corner point estimates for each of the nc components in the mixture. The mole fraction of the corner point component in the trial phase is set to 0.9. Li and Firoozabadi22 noted that stability of a two-phase mixture could be assessed through stability testing of the higher molecular weight phase. The use of several initial K-value estimates is necessary for the detection of phase instability in nonpolar hydrocarbon fluids, particularly in the near-critical region.10,22 In the case of hydrocarbon−water mixtures, fewer initial K-value estimates may suffice. Some researchers have suggested using just three initial K-value estimates, corresponding to the heaviest component, the ideal gas law, and a correlation for water solubility in the hydrocarbon phase.23 However, other authors have maintained the approach in eq 2 proposed by Li and Firoozabadi.22 Mahabadian et al.24 used Nc + 4 initial K-value estimates for stability test calculations in multiphase hydrate mixtures. Most recently, Pang and Li25 used the same Nc + 4 initial estimates in addition to a near-pure water estimate for the trial phase composition. A recently developed algorithm for hydrocarbon−water mixtures used nested two-phase flash calculations and stability tests, including a hydrocarbon liquid−vapor stability test assuming no water was present in the system.26 However, this method used Henry’s Law and was focused on CO2-rich systems at moderate temperatures. To the best of our knowledge, there is no dedicated protocol for stability testing of hydrocarbon−water mixtures in a thermal setting. Phase-behavior modeling can be particularly arduous in reservoir simulation, especially for three-phase systems in models with a large number of grid blocks. One potential solution is the use of reduced-order models for phaseequilibrium calculations. Reduced-order models (herein referred to as reduced-variables methods) involve recasting an underlying problem with a smaller set of independent variables. Reduced variables were first used for phase-behavior computations in the early 1980s.27 Phase-split calculations were first solved with reduced variables using an assumption of all zero-BIPs.28 A generalized method was then introduced by Hendriks29 and later extended by Hendriks and Van Bergen.30 Additional developments in the use of reduced variables have been made by several researchers since, including the use of new sets of parameters for stability-analysis and phase-split calculations.31−36 Nichita et al.37 were the first to extend reduced-variables methods to multiphase equilibrium. Despite the ubiquity of reduced-variables formulations, there is only one implementation of a two-phase reduced-variables method in a commercial reservoir simulator.38,39 Several recent studies have compared the performance of reduced variables with conventional methods when used for flash calculations.40−43 Michelsen et al.41 observe that performance is highly implementation dependent, and an objective comparison is further complicated by the fact that some reduced-variables formulations solve the problem using a reduced basis, thereby solving what is a different problem from the outset. Although there have been several comparison studies with conventional methods, the application of reduced formulations has been restricted almost exclusively to hydrocarbon mixtures. One notable exception is the use of reduced variables for four-phase flash calculations accounting for CO2 solubility in water.44 However, the application was restricted to an isothermal setting involving light hydrocarbon components and a highly water-soluble injectant. In addition, the reduced

method involved modification of the foundational mixing rules used in the equation of state. Reduced-variables methods have not been investigated for phase-equilibrium computations of hydrocarbon−water mixtures in which trace components appear or for high-temperature applications where water appears in the vapor phase. The focus of this work is a general purpose methodology for phase-equilibrium calculations involving hydrocarbon−water mixtures in thermal reservoir simulation. Our research consists of two contributions: I. Physics-based phase-stability testing. II. Development and implementation of a reduced-variables method for three-phase equilibrium calculations. This paper is organized as follows. First, we review the key theory on phase-stability testing and phase-split calculations, including the numerical solution of these problems. Second, we present our solution strategy. We use the distinctive behavior of hydrocarbon−water systems that allows us to develop a robust stability analysis procedure and that makes phaseequilibrium calculations amenable to the reduced method. We extend the method of Pan and Firoozabadi32 to three phases. Next, we examine the performance of our approach over a wide range of pressures and temperatures using several case studies of hydrocarbon−water mixtures. Finally, we examine the reasons for the superior performance of the reducedvariables approach for hydrocarbon−water mixtures. Based on the results of our research, we advocate for the use of the reduced-variables approach to rigorous multiphase-equilibrium calculations for hydrocarbon−water mixtures in thermal reservoir simulation.

■

BACKGROUND Stability Analysis. The concept of tangent plane distance was introduced by Baker et al.45 and formally expressed by Michelsen.17 Instability of a test phase, z, is indicated by a negative tangent plane distance: Nc

TPD(yi ) =

∑ yi (ln ϕi(̂ yi ) + ln yi − ln ϕi(̂ zi) − ln zi) i=1

(3)

where ϕ̂ i refers to the fugacity coefficient of component i. Fugacity coefficient is a measure of nonideality and is defined as the ratio of fugacity to its value at the ideal state ϕ̂ i = f i/ (yiP). Note that y is the vector of mole fractions of a trial phase emanating from the test phase. Stability analysis rests on the premise that if all stationary points of the TPD function are located, the global minimum will be located. In practice, this is very difficult because the Gibbs energy surface is inherently nonconvex and practical solution methods converge to local stationary points. Michelsen17 introduced a new set of primary variables, Yi = yi × exp(μi(z) − μi(y)). Formally, the new independent variables Yi are mole numbers, with the corresponding mole fractions given by yi = Yi/∑jYj. The stationary condition can be written as ln(Yi ) + ln ϕi(̂ Yi ) − ln ϕi(̂ zi) − ln zi = 0

(4)

Stationary points correspond to solutions of eq 4. Stability is verified, provided that at all stationary points TPD ≥ 0, which N is equivalently, ∑i cYi ≤ 1. Michelsen17 suggested an equivalent, unconstrained version of the minimization problem as follows: C


Article

Industrial & Engineering Chemistry Research Nc

TPD*(Yi ) = 1 +

Nc

∑ Yi(ln ϕi(̂ Yi) + ln Yi − ln ϕi(̂ zi) − ln zi − 1)

∑

i=1

i

(5)

V

fi ̂

W

− fi ̂

L

− fi ̂

=0

(6)

=0

(7)

L

∑ i

(8)

ln K iW = ln ϕ(̂ xi) − ln ϕ(̂ wi)

(9)

ϕ(̂ xi) ϕ(̂ y )

(10)

ϕ(̂ xi) ϕ(̂ w )

(11)

i

K iW =

i

zi(K iW − 1) 1 + V (K iV − 1) + W (K iW − 1)

=0 (12)

=0 (13)

Compositions of the remaining phases are subsequently found through connection to the reference phase by way of the relevant equilibrium ratios: yi = K iVxi

(15)

wi = K iW xi

(16)

The choice of independent variables influences the solution strategy. Popular choices are component mole numbers, equilibrium ratios, and the logarithm of equilibrium ratios. Regardless of the independent variables, flash calculations are typically initiated via solution of the nonlinear equations using successive substitution iteration. A Newton method may be used if sufficiently close to the solution. When the K-values or the natural logarithms of the K-values are used, the nonlinear equations are solved via Newton’s method. In the case of component mole numbers, the nonlinear equations may be solved via the Newton method or the problem may be formulated for minimization of Gibbs free energy. In this research we use the Newton method, which is efficient, although unlike the minimization approach there is no guarantee of convergence. Michelsen51 noted that the most troublesome elements of any isothermal flash algorithm are the initial estimate of Kvalues and attaining convergence in the near-critical region. In the context of reservoir simulation, good initial estimates of phase-equilibrium ratios may be obtained from previous time steps or the previous iteration. With small CFL numbers (for instance, in IMPES implementations) or far from displacement fronts, initial estimates from prior time steps may be an excellent starting point for the phase-split calculation. However, phase changes may occur very frequently in reservoir simulation, such as in the vicinity of displacement fronts or when large time steps are taken in fully implicit implementations. In thermal reservoir simulation these phase changes are also driven by temperature gradients. Convergence in the critical region is an important problem to address in the case of hydrocarbon mixtures, particularly in the study of phase equilibria of miscible displacements. In the case of thermal recovery processes, the operating conditions are generally far from criticality. However, bicritical points may still be encountered, whereby two hydrocarbon phases approach criticality while the water component exists in a near-pure aqueous phase.

where yi, xi, and wi refer to mole fractions of component i in the vapor, oleic, and aqueous phases, respectively. Rearranging these equations, yields the definition of the equilibrium ratios (K-values) in terms of fugacity coefficients: K iV =

− 1) + W (K iW − 1)

where zi is the mole fraction of component i in the feed composition, and V and W are the mole fractions of the vapor and aqueous phases, respectively. Here we take the oleic (liquid hydrocarbon) phase as the reference phase. Equations 12 and 13 yield phase amounts upon convergence. Component mole fractions of the reference oleic phase are then computed as a function of the phase-equilibrium ratios and phase amounts: zi xi = 1 + V (K iV − 1) + W (K iW − 1) (14)

where f ̂ refers to fugacity, and the superscripts V, L, and W refer to the vapor, oleic and aqueous phases, respectively. The equality of component fugacities is often expressed using the logarithm of the equilibrium ratios and phase fugacity coefficients:49 ln K iV = ln ϕ(̂ xi) − ln ϕ(̂ yi )

1+

Nc

Here, Yi > 0 is the only restriction. The same stationary points with the same sign are identified with TPD and TPD*. Stability analysis is always performed as a test for instability of a single phase, regardless of the number of phases present. Multiphase-stability analysis is performed through a stepwise approach of a single-phase stability test, followed by two-phase flash, followed by single-phase stability testing of one of the two equilibrium phases. The two-phase stability test is thermodynamically equivalent to single-phase stability testing. In reservoir simulation, phase-equilibrium computations are performed across many thousands of grid blocks at each time step and (in some formulations) at each iteration. In compositional simulation, stability analysis becomes a significant computational burden. Several authors have proposed strategies to bypass stability analysis.46−48 This is not the focus of this paper. Phase-Split Calculations. For a given number of phases, the phase-split calculation solves both for phase compositions and phase amounts. At convergence, the phase-split calculation yields the minimum Gibbs free energy of the system. This is achieved by satisfying the constraints of equality of component fugacities across phases and material balance equations. For brevity, we exclude the development of the two phase-split problem here and focus on the three-phase system. For a threephase hydrocarbon−water system, the isofugacity constraints are expressed as

fi ̂

zi(K iV − 1) V (K iV

The fugacity coefficients are computed as a function of pressure, temperature, and composition, using the cubic EOS. In addition to the isofugacity constraints, the component material balance constraints must be satisfied. Mole balance constraints are easily expressed in terms of the Rachford−Rice objective functions:50 D


Article

Industrial & Engineering Chemistry Research

■

FORMULATION AND SOLUTION STRATEGY Phase-equilibrium calculations for three-phase (VLL) hydrocarbon systems have been studied for over 30 years.10,17,22,51 These phase-equilibrium calculations involve crude oil and CO2 or other miscible gases. The high degree of miscibility leads to complex phase behavior. A comprehensive set of Nc + 4 initial K-value estimates has been used for stability-testing calculations in these mixtures.22,52 On the other hand, because water and liquid hydrocarbon phases are essentially immiscible, the phase behavior of the three-phase VLW system is far more predictable than that of highly miscible VLL systems. For twophase hydrocarbon vapor−liquid systems, the initial estimates Wilson and inverse Wilson are sufficient for stability testing.10 For water−hydocarbon systems an additional immiscible aqueous phase may exist, to which we assign an initial estimate KH2O. Three estimatesWilson, inverse Wilson, and KH2O are sufficient for phase-stability tests in VLW systems. The use of Nc + 4 calculations for both single and two-phase stability testing in VLW systems is impractical. Three-phase VLL and VLW systems are of very different miscibility. We exclude the Nc + 4 initial K-value estimates in phase-stability testing calculations involving VLW systems. We only use the three initial estimates: Wilson, inverse Wilson, and KH2O. Objectives. We now present a solution strategy for phaseequilibrium calculations for hydrocarbon−water mixtures. We perform stability-analysis and phase-split calculations in sequence. This is the stepwise procedure originally suggested by Michelsen.51 Our objective here is twofold:

Figure 1. Partitioning of the pressure−temperature domain using the steam saturation curve. The aqueous phase does not exist below the steam saturation pressure. Above the steam saturation pressure an aqueous phase may be present, and three sets of initial K-values are 2O , 1/KWilson , KH }. used for stability analysis, {KWilson i i i

where τ = (1 − T/Tc). Equation 18 provides a guess that is very close to steam saturation pressure computed with the PR EOS, and the cost of the iterative solution is very small. The fugacity of pure water is a function of pressure and temperature and is calculated using the PR-EOS. Once the temperature is specified, the steam saturation pressure is calculated by solving the nonlinear equation

I. Development of a general purpose protocol, accounting for all components in all phases (including trace quantities).

f wV − f wW = 0

f Vw

where and refer to fugacity of water in the vapor and aqueous phases, respectively. We use the secant method to solve for pressure. The initial estimate for steam saturation pressure is found using eq 18. Since the initial value is very close to the solution, less than five iterations are typically required to converge to a residual of 10−6. The correlation in eq 18 is accurate and widely used in industry. It serves as a good initial guess for steam saturation pressure, prior to solution of eq 19 with the PR EOS. If pressure is below the pressure of saturated steam, only the Wilson and inverse-Wilson correlations are used to initiate stability analysis. If pressure is greater than the saturation 2O pressure of steam, KH is used as an additional initial guess. i The physical reasoning underpinning this strategy is that the aqueous phase may exist only at pressures above the steam saturation pressure. If phase instability is detected, two-phase flash calculations are performed. A simple correlation for water solubility in oil54 is used to initiate the two-phase split calculation when instability of the mixture is detected above steam saturation pressure:

II. Consistent resolution of phase boundaries and convergence of phase-split calculations to tight tolerance for reservoir simulation. Stability Analysis. We propose a scheme that employs a maximum of three sets of initial guesses for single-phase stability testing: Wilson, inverse-Wilson, and near-pure H2O: {KIinitial} = {K iWilson , 1/K iWilson , K iH 2O} HO

(19)

fW w

(17) HO

where Ki 2 = 0.999/ztest for i = H2O, and Kj≠i2 = 0.001/(Nc − i test 1)/zj . Our rationale is that the appearance of a trial phase predominantly comprised of hydrocarbons will be detected by stability tests initiated with the Wilson correlation or its inverse. On the other hand, the evolution of an aqueous or near-pure steam trial phase will be detected by stability tests HO initiated using Ki 2 . The insight here is that the saturated-steam curve can be used to guide the selection of trial phases. This approach is sensitive to the distinct behaviors of hydrocarbon−water mixtures. We begin by partitioning the pressure−temperature domain, as illustrated in Figure 1. In the interests of thermodynamic consistency, we calculate the steam saturation pressure Psat using the PR EOS.12 This requires iteration, for which we use the secant method. We use a correlation53 to provide an initial guess for Psat:

ln(x HL2O) = −21.2632 + 5.9473 × 10−2T − 4.0785 × 10−5T 2 (20)

xLH2O

where is the mole fraction of water in the oleic phase and L T is the temperature in Kelvin. To determine KH2O from xH 2O 23 we follow the approach outlined in Patterson et al. Assuming the aqueous phase is essentially water, ln(xW H2O) = 0. Thus, L W L ̂ ̂ ln(ϕH2O) = ln(ϕH2O) − ln(xH2O). The K-value is simply calculated as the ratio of fugacity coefficients. Consistent with Raoult’s law, the bubble point pressure of a reservoir aqueous phase containing dissolved hydrocarbons

ÅÄÅ T P sat = Pc × 10 × expÅÅÅÅ c ( −7.85823 + 1.83991τ 1.5 − 11.7811τ 3 ÅÅÇ T ÉÑ Ñ + 22.6705τ 3.5 − 15.9393τ 4 + 1.77516τ 7.5)ÑÑÑÑ ÑÑÖ (18)

E


Article

Industrial & Engineering Chemistry Research and nonhydrocarbons (CO2, H2S, N2) is higher than that of an aqueous phase consisting of pure water. When the system pressure P is less than the saturated pressure (Psat) of a pure water phase at a given temperature, there is no aqueous phase present in the system. For cases in which P > Psat, a second two-phase split calculation is performed using the set of K-values obtained from the stability test corresponding to the stationary point with the lowest negative TPD. The stability of the two-phase mixture is then assessed. The two-phase split result yielding the lowest Gibbs energy is first selected for stability analysis. The phase state of the two-phase system guides the choice of the next set of K-values. If the aqueous phase is present after the two-phase split calculation, we use the Wilson correlation and its inverse to detect instability of the hydrocarbon phase. If no aqueous phase is present, the higher molecular weight phase is tested for stability using KH2O as the starting set of K-values. A flowchart of the approach to stability analysis is included in Appendix A, and a stepwise description is provided below. • Step 1: Specify P, T, and z i , in addition to thermodynamic parameters for each feed component including molecular weight, Tc, Pc, ω, and binary interaction terms. • Step 2: Calculate the steam saturation pressure, Psat, at the specified temperature. • Step 3: If P < Psat, the aqueous phase cannot be present. Perform stability analysis using the Wilson correlation and its inverse and proceed to step 4. If P > Psat, the aqueous phase may be present. Perform stability analysis using three sets of K-values {KWilson, 1/KWilson, KH2O} and proceed to step 5. • Step 4: (a) If a nontrivial negative TPD is found, perform a two-phase flash. Initiate flash with the set of K-values from stability testing that yields the minimum negative TPD. Stop upon convergence of two-phase flash. (b) If no negative TPD is found, the system is single phase. Stop. • Step 5: (a) If a nontrivial negative TPD is found, perform two separate two-phase flash calculations. For the first flash, use the correlation in eq 20 to determine the fraction of H2O in the oleic phase and estimate KH2O. For the second flash, use the set of K-values from the stability test corresponding to the minimum negative TPD. Proceed to step 6. (b) If no negative TPD is found, the system is single phase. Stop. • Step 6: If the two-phase flash results are identical and more than one negative TPD was detected in stability analysis, perform an additional two-phase flash using the set of K-values yielding the next lowest TPD. Select the two-phase flash result with the lowest Gibbs energy and proceed to step 7. • Step 7: Test for stability of the two-phase state. (a) If the aqueous phase is present, test the stability of the hydrocarbon phase using the Wilson correlation and its inverse. (b) If the aqueous phase is not present, test for the appearance of an aqueous phase using the phase with higher molecular weight. • Step 8: (a) If instability is detected, perform three-phase flash. To initialize the flash, use the K-values from twophase flash and the K-values from the stability test in step 7. (b) If instability is not detected and an additional unique two-phase flash result is available, select this two-

phase flash result and go back to step 7. Otherwise, the system is two-phase. The stability-analysis strategy has two direct benefits. First, the number of trial phases required for stability testing is predetermined based on the system pressure and temperature. Second, the reliability of the phase-stability tests and the proceeding phase-split calculations is greatly improved. Location of TPD function (eq 5) stationary points is the overarching difficulty in stability analysis, and so the focus of our algorithm is strategic rather than numeric. We use the same calculation procedure with reduced variables as outlined in Firoozabadi and Pan.31 The stability test calculation is far simpler than the solution of the phase-split problem because a stability test is not subject to material balance constraints and only involves a single variable phase composition. We now proceed to describe the numerical solution of the three-phasesplit calculation. Three-Phase-Split Using a Reduced Variables Method. The most pressing difficulties in flash calculations are obtaining initial estimates and convergence in the critical region.51,55 To be clear, we use the K-values from stability analysis to address the former challenge. Convergence in the critical region has not been addressed specifically for hydrocarbon−water mixtures to the best of our knowledge. Hydrocarbon−water mixtures are highly immiscible at low temperatures. However, in thermal EOR processes this is not the case. Successive substitution iteration (SSI) and the Newton method are typically used for phase-split calculations. Newton methods require the solution to be sufficiently close to the true solution for convergence. For this reason, SSI is used before switching to the Newton method once the residual norm of the fugacity equations (eqs 8 and 9) is sufficiently small. Michelsen51 noted that the strongly nonideal behavior of systems involving liquid hydrocarbon and liquid water phases leads to convergence of successive substitution within a few iterations. Aqueous−liquid phase coexistence occurs in the vast majority of phase-split calculations at low temperatures. In a thermal setting, multiple phase states are possible and critical regions can be encountered. Successive substitution is ill-suited to the solution of phase-split calculations in the critical region, as K-values approach unity. Convergence of successive substitution becomes prohibitively slow as the critical point is approached. Either the Newton method or a quasi-Newton method is required to reach the solution in a reasonable number of iterations. We extend the method of Pan and Firoozabadi32 to three phases. We use the PR EOS,12 which for a pure component is written P=

ai(T ) RT − v − bi v(v + bi) + bi(v − bi)

(21)

where R is the gas constant, b denotes the covolume, and a refers to the energetic parameter. The crux of reduced variables methods lies in the decomposition of the energetic parameter, which for a cubic EOS such as the PR EOS is given by Nc

a=

Nc

∑ ∑ xixj i=1 j=1

F

aiaj (1 − σij) (22) DOI: 10.1021/acs.iecr.9b00695 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research where the binary interaction parameter (BIP) matrix is β = (1 − σij). Given that the BIP matrix is a square, symmetric matrix, the eigenvalue decomposition yields orthogonal basis vectors.30 A = QΛQ−1 = QΛQT

• Step 4: Calculate the independent phase amounts V and W, by solving the three-phase Rachford−Rice equations (eqs 12 and 13). Update phase compositions xi, yi, wi via eqs 14, 15, and 16. • Step 5: Check the fugacity residual ∥fL̂ − fV̂ ∥ + ∥fL̂ − ̂ ∥. If the residual is smaller than the tolerance in the fW convergence criterion, stop. If the residual is smaller than the switching criterion tolerance, proceed to the Newton procedure. Otherwise, update the reduced variables {QV, QW, bV, bW, V, W} and go back to step 2. The Newton procedure is as follows: • Step 1: Initialize the reduced variables {QV, QW, bV, bW, V, W} • Step 2: Calculate {QL, bL, L}. • Step 3: Calculate fugacity coefficients, ϕ̂ Li , ϕ̂ Vi , ϕ̂ W i and ̂L ̂W the equilibrium ratios KVi = ϕ̂ Li /ϕ̂ Vi and KW i = ϕi /ϕi . • Step 4: Calculate phase compositions xi, yi, wi using eqs 14, 15, and 16. • Step 5: Check the fugacity residual ∥fL̂ − fV̂ ∥ + ∥fL̂ − ̂ ∥. If the residual is smaller than the tolerance in the fW convergence criterion, stop. Otherwise, proceed to step 6. • Step 6: If convergence is not achieved, construct and solve the linear system

(23)

where Q is an orthogonal matrix, with column vectors corresponding to the eigenvectors of β: Q = {q′ (1) , q′ (2) , ..., q′ (Nc)}

(24)

and Λ is a diagonal matrix with entries corresponding to the eigenvalues of β: ÉÑ ÄÅ ÑÑ ÅÅ λ1 ÑÑ ÅÅ ÑÑÑ ÅÅÅ ∏ Λ = ÅÅ ÑÑ ÑÑ ÅÅ Ñ ÅÅ λ ÑÑÖ ÅÅÇ Nc Ñ (25) Each eigenvector may be written q′(i) = (qi1′ , qi2′ , ..., q′iNc)T. By extension Nc

β=

∑ λkqki′ qkj′

(26)

k=1

Combining eq 22 and 26 yields Nc

a=

Nc

J·ΔX = −e

Nc

∑ λk ∑ xiqki′ ai

∑ xjqkj′ aj

k=1

j=1

i=1

• Step 7: Update the reduced variables and go back to step 2. The solution of the three-phase Rachford Rice equations (eqs 12 and 13) is referred to as a constant K-value flash.56 The constant K-value flash is an integral component of SSI in the solution of the phase split. In this research we use a trustregion optimization solver to minimize a convex objective function, similar to the line-search solution developed by Okuno et al.57 There exists a trade off in choosing the tolerance in the switching criterion. Switching too early from successive substitution to Newton’s method can result in failure of the Newton method. Switching too late can result in excessive successive substitution iterations. We choose a default switching criterion of 0.02. If the Newton method fails, we divide the switching criterion tolerance by a factor of 10 and resume successive substitution iteration. Failure of Newton’s method is indicated by a negative phase amount. The derivation and implementation of analytical derivatives is challenging and error prone. To verify the integrity of our implementation, we have compared the partial derivatives of the Jacobian matrix with those generated using the Stanford Automatically Differentiable Expression Templates Library (ADETL).58,59 The full development of the nonlinear equations and their derivatives is provided in Appendix B. The equilibrium state corresponds to the global minimum Gibbs free energy. Phase-split calculations can converge to a local minimum Gibbs free energy. In our implementation we compute the normalized Gibbs free energy G̅ after two- and three-phase split calculations. For a two-phase mixture, the normalized Gibbs free energy is

(27)

Now, letting qki = qki′ ai and defining Qk = ∑i N=c 1xiqki, we can rewrite a in reduced form: Nc

a=

∑ λkQ k2

(28)

k=1

In a cubic EOS the compressibility factor and fugacity coefficient can be written as a function of the energy and covolume parameters: Z = Z(Q 1 , Q 2 , ..., Q N , b)

(29)

ϕî = ϕi(̂ Q 1 , Q 2 , ..., Q N , b)

(30)

c

c

However, in many cases sparsity in the BIP matrix yields eigenvalues that are close to zero or insignificant, and both Z and ϕ̂ i become functions of a reduced set of parameters: Z = Z(Q 1 , Q 2 , ..., Q m , b)

(31)

ϕî = ϕi(̂ Q 1 , Q 2 , ..., Q m , b)

(32)

(33)

where m < Nc. In a three-phase setting, the set of reduced variables is {QV, QW, bV, bW, V, W}. The reduced variables method acts as a low-pass filter. This ensures a full-rank linear system, and the Newton method can be used. A stepwise description of the three-phase flash procedure is provided here for clarity. A description of the two-phase flash is available in Pan and Firoozabadi.32 The successive substitution iteration is as follows: • Step 1: Initialize the reduced variables {QV, QW, bV, bW, V, W}. • Step 2: Calculate {QL, bL, L}. • Step 3: Calculate fugacity coefficients, ϕ̂ Li , ϕ̂ Vi , ϕ̂ W i and ̂L ̂W the equilibrium ratios KVi = ϕ̂ Li /ϕ̂ Vi and KW i = ϕi /ϕi .

Nc

G̅ = G /RT =

∑ (niL ln fî

L

V

+ niV ln fi ̂ )

i

(34)

By extension, for a three-phase system G


Article


G̅ =

∑ (niL ln fî

L

+ niV ln fi ̂

V

using the reduced formulation proceed rapidly, in spite of the relatively large number of pseudocomponents. This case is an example of the difficulty in identifying a narrow three-phase state. Figure 2 shows the narrow threephase region in pressure−temperature space. We use the phase equilibrium calculation at 85 bar and 571 K as an example. This point is above the steam saturation curve and stability analysis yields a negative TPD, indicating instability of the mixture. Table 3 shows the results of a two-phase flash calculation initiated with K-values predicted using eq 20 and assuming an almost pure aqueous phase in combination with a liquid hydrocarbon phase. Table 4 shows the results of an alternative two-phase split calculation, initiated using K-values predicted by stability analysis. The normalized Gibbs free energy is calculated using eq 34 for two-phase systems. The two-phase system in Table 4 has a lower normalized Gibbs free energy than the two-phase system in Table 3. However, stability testing of the two-phase mixture in Table 4 does not indicate the appearance of a third phase. Stability testing of the hydrocarbon phase in Table 3 does yield a negative TPD, and the results of the ensuing three-phase split calculation are displayed in Table 5. For three-phase systems we compute normalized Gibbs free energy using eq 35. Case 2: Four-Component Heavy Oil Plus Water. Case 2 is taken from the dissertation of Brantferger.62 This case is used to demonstrate the efficacy of the stability testing protocol in narrow regions of immiscibility. Table 6 lists properties of the four-component heavy oil and water. Binary interaction parameters are listed in Table 7. The spectral decomposition of the BIP matrix yields only two nonzero eigenvalues. Figure 3 shows that as temperature is increased the equilibrium state changes from two-phase liquid−aqueous (LW), to three-phase vapor−liquid-aqueous (VLW), and finally to two-phase vapor−liquid (VL). Delineation of the three-phase region is difficult, and ordinary stability testing protocols are inadequate. Single-phase stability analysis yields two negative TPD results. Given that the mixture is unstable, we perform two-phase flash using the K-values predicted via eq 20 and those generated from the stability test yielding the lowest negative TPD. Our approach to stability testing involves multiple two-phase flash calculations. This is important for correct identification of the phase state. The generalized strategy for multiphase stability testing and flash described by Li and Firoozabadi22 consists of Nc + 4 initial K-value estimates (see eq 2). In this approach, failure of a three-phase flash calculation is used as an indication to perform an additional two-phase flash calculation with the set of K-values using the next lowest negative TPD from single-phase stability testing. The higher molecular weight phase from the second two-phase flash result is then tested for instability, and a three-phase flash is performed if necessary. While this method was demonstrated for several complex hydrocarbon fluid mixtures and hydrate-containing mixtures,24 it has been shown to fail10 and is inadequate for this case. Our approach circumvents difficulties in phase identification by seeking a unique two-phase flash result. The results of our approach to stability testing are displayed in Figure 4. Figure 5 displays the system depicted in Figure 4 in pressure−enthalpy space. Here, a very small three-phase region in the PT-phase diagram corresponds to a large region in PH space. Meticulous resolution of phase boundaries and

W

+ niW ln fi ̂ )

i

(35) 32

Although we extend the method of Pan and Firoozabadi to three phases for the first time, application of reduced variables methods beyond two-phase flash has been previously established. Several authors have presented reduced formulations for three- and four-phase equilibrium.37,44,56,60 However, our analysis in the proceeding section examines the use of reduced variables for hydrocarbon−water mixtures in a thermal setting. This is the first application of reduced variables to thermal compositional systems. In addition, the research described herein is the first to examine the applicability of a reduced variables method to systems characterized by the presence of trace components in certain phases.

■

RESULTS AND ANALYSIS In this section we examine the performance of both our stability analysis strategy and the reduced-variables formulation for two- and three-phase-split calculations. We begin by examining the performance of the stability-analysis algorithm in delineating phase boundaries for hydrocarbon fluids in the presence of water at elevated temperatures and pressures. We tested our method across a large parameter space. Compositional variation was achieved by varying the fraction of H2O combined with the hydrocarbon mixtures. We altered the fraction of water in the mixture from 1% to 99%, with consistent resolution of phase boundaries observed across the compositional spectrum. Here we only show the most material test cases to demonstrate the efficacy of our approach. Our analysis shows that the stability-testing strategy reliably identifies phase instability and leads to the correct result. In addition, the reduced-variables method performs consistently in solving the phase-split problem. Case 1: Seven-Component Heavy Crude Plus Water. Case 1 is taken from Lapene et al.61 The oil is a heavy Venezuelan crude from a thermal recovery project. The composition is provided in Table 1, with binary interaction Table 1. Properties of Seven-Component Heavy Oil Plus Water Mixture, Taken from Lapene et al.61,a component

mole fraction

MW (g/mol)

Tc (K)

Pc (bar)

ω

H2O C2−C11 C12−C16 C17−C21 C22−C27 C28−C35 C36−C49 C50+

0.99 0.0001807 0.0018070 0.0015660 0.0013850 0.0011450 0.0010240 0.0028923

18.015 143.66 193.82 263.40 336.29 430.48 573.05 1033.96

647.37 635.64 701.24 772.05 826.30 879.55 936.97 1260.0

221.20 24.11 19.25 15.10 12.29 9.94 7.79 6.00

0.344 0.4645 0.6087 0.788 0.9467 1.1042 1.273 1.65

a

The hydrocarbon components are characteristic of a heavy crude.

parameters listed in Table 2. This is a complex hydrocarbon fluid that has been lumped into seven pseudocomponents. A very small fraction (less than 2%) of the crude is comprised of intermediate hydrocarbon components, represented by a single pseudocomponent (C2−C11). Characterization with a practical number of pseudocomponents is challenging for heavy oils. The spectral decomposition of the BIP matrix yields only two significant eigenvalues. As a result, the phase-split calculations H


Article


Table 2. Binary Interaction Parameters for Mixture of Seven-Component Heavy Crude with Water, Taken from Lapene et al.61 BIP

H2O

C2−C11

C12−C16

C17−C21

C22−C27

C28−C35

C36−C49

C50+


0 0.0952 −0.48 −0.48 0.45 0.53 0.50 0.50

0 0 0 0.13 0.135 0.1277 0.1

0 0 0.05 0.08 0.1002 0.1

0 0 0 0.09281 0.130663

0 0 0 0.006

0 0 0.006

0 0

0

Table 5. Three-Phase Phase Flash Result Obtained after Testing the Stability of the Two-Phase System in Table 3a component

phase 1

phase 2

phase 3


0.995912 0.000214512 0.00194683 0.00120972 0.000547622 0.000144005 2.50245e−005 1.05125e−009

0.237088 0.00133969 0.030944 0.0690683 0.107303 0.116028 0.112949 0.32528

1 2.39218e−010 2.40812e−011 1.86274e−014 4.44569e−018 4.16862e−023 9.59619e−031 6.03852e−051

a

Stability testing of the two-phase system in Table 4 did not indicate phase instability. G̅ = −0.360728432276.

Table 6. Properties of the Five-Component Mixture in Case 262,a

Figure 2. Mixture of 99% water and 1% oil, as described in Table 1. Pressure is sampled in increments of 0.1 bar, and temperature is sampled in increments of 0.05 K. The three-phase region is very narrow in pressure−temperature space.

component

mole fraction

MW (g/mol)

Tc (K)

Pc (bar)

ω

H2O C8 C13 C24 C61+

0.50 0.2227 0.1402 0.1016 0.0355

18.015 116.00 183.00 337.00 858.00

647.37 575.78 698.00 821.30 1010.056

221.20 34.82 23.37 12.07 7.79

0.344 0.400 0.840 1.070 1.330

Table 3. Two-Phase Flash Results for Heavy Oil Plus Water Mixture in Case 1 at 85 bar and 571 Ka component

phase 1

phase 2


1 2.41968e−09 1.05679e−10 3.18043e−14 4.32222e−18 3.09172e−23 6.48657e−31 4.21787e−51

0.244793 0.0136464 0.136466 0.118265 0.104596 0.0864712 0.0773332 0.218429

The fluid consists of water in combination with a heavy oil. The oil is devoid of light hydrocarbon components. a

Table 7. Binary Interaction Parameters for Five-Component Mixture in Case 262,a

G̅ = −0.3586201422.

a

Table 4. Two-Phase Flash Results for the Mixture in Case 1a

BIP

H2O

C8

C13

C24

C61+

H2O C8 C13 C24 C61+

0 0.5 0.5 0.5 0.5

0 0 0 0

0 0 0

0 0

0

a

component

phase 1

phase 2


0.996479 0.000172941 0.00160261 0.00105735 0.000517883 0.000144512 2.56697e−005 1.06639e−009

0.236495 0.00108307 0.025577 0.0607212 0.10223 0.117501 0.117129 0.339264

The spectral decomposition of the BIP matrix results in only two non-zero eigenvalues, owing to the large number of zero interaction coefficients.

reservoir simulation. Figure 6 clearly shows the sharp change in enthalpy at the phase boundary. Case 3: Three-Component Oil Plus Water. Case 3 is taken from Varavei.63 Table 9 lists properties of the threecomponent oil in combination with water. Binary interaction parameters are listed in Table 10. This case study introduces an entirely different fluid system from the heavy crude oils of the first two case studies. The full spectrum of two-phase states (LW, VW, VL) is present in Figure 7. The vapor−liquid region (in yellow) extends above the steam saturation curve, in spite of the high concentration of water in the system. Correct

a

Results are from the second two-phase split calculation, performed at 85 bar and 571 K. G̅ = −0.3606695361.

quantification of phase fractions is important for the fidelity of the coupling to the energy and material balance equations in I


Article

Industrial & Engineering Chemistry Research Table 8. Heat Capacity Coefficient Data for the Fluid in Case 262 component H2O C8 C13 C24 C61+

C0P1 J/(mol K) 32.200000 −1.23e+01 −5.080000 −5.690000 0.1230000

C0P2 J/(mol K)

C0P3 J/(mol K)

C0P4 J/(mol K)

−05

−3.596 × 10−09 0.0 0.0 0.0 0.0

1.055 × 10 −2.52e−04 −4.14e−04 −7.64e−04 −1.95e−03

0.001924 6.65e−01 9.97e−01 1.840000 4.750000

Figure 5. Mapping of Figure 4 from the pressure−temperature domain to pressure−enthalpy space. An ostensibly narrow three-phase region maps to a large region in pressure−enthalpy space.

Figure 3. PT-phase diagram for the 50% water/50% oil mixture in Case 2. Component properties are provided in Table 6. Pressure is sampled in increments of 0.1 bar, and temperature is sampled in increments of 0.1 K. The mixture is characterized by a narrow threephase region.

Figure 6. Enthalpy contours corresponding to the PT-phase diagram in Figure 3. At the phase boundary there is an abrupt change in enthalpy.

Figure 4. Close up of PT-phase diagram shown in Figure 3, showing consistent resolution of the phase state. Pressure is sampled at 0.1 bar intervals from 2 bar through 20 bar. Temperature sampled at 0.1 K intervals from 380 to 480 K.

Table 9. Component Properties for Four-Component Hydrocarbon−Water System in Case 363,a

identification of the phase state is not straightforward because variations in pressure, temperature, and composition cause phase transitions across several two-phase and three-phase states (LW, LV, VW, VLW). Resolving the two-phase vapor− aqueous region (VW) is particularly challenging, as the water component is present in high concentrations in the vapor phase. The oft-used approach of initiating flash computations assuming a two-phase LW system is inappropriate for this system, particularly at temperatures close to steam saturation temperature. The three-phase region in Figure 7 is narrow, as observed in the first two case studies. We emphasize that the uncertainty regarding phase state is at pressures above the

component

mole fraction

MW (g/mol)

Tc (K)

Pc (bar)

ω

H2O C6 C10 C15

0.950 0.005 0.015 0.030

18.015 86.178 134.00 206.00

647.3 507.5 622.1 718.6

220.4732 32.88847 25.34013 18.49090

0.344000 0.275040 0.443774 0.651235

a

Compared to Cases 1 and 2 the hydrocarbon components are relatively volatile.

steam saturation curve. At pressures below the steam saturation curve, the Wilson correlation and its inverse are sufficient to detect the instability of the mixture. A vapor− liquid state is the only possible two-phase state at pressures J


Article


Okuno3 are listed in Table 12. For this fluid there are only two significant eigenvalues, as the BIP matrix consists primarily of

Table 10. Binary Interaction Parameters for FourComponent Mixture in Case 363 BIP

H2O

C6

C10

C15

H2O C6 C10 C15

0 0.48 0.48 0.48

0 0.002866 0.010970

0 0.002657

0

Table 12. Binary Interaction Parameters for the Water− Pseudocomponent Oil Mixture in Table 113,79 BIP

H2O

PC1

PC2

PC3

PC4

H2O PC1 PC2 PC3 PC4

0 0.71918 0.45996 0.26773 0.24166

0 0 0 0

0 0 0

0 0

0

zeros. Figure 8 shows the phase diagram for a mixture of 50% water and 50% oil. The three-phase region is very large for this system, particularly in comparison to the previous case studies.

Figure 7. PT-phase diagram for the fluid in Case 3. The system is a mixture of 95% water and 5% oil. Component properties are provided in Table 9. Multiple two-phase states (VW, VL) border the narrow three-phase region. Pressure is sampled in increments of 0.25 bar, and temperature is sampled in increments of 0.25 K.

below steam saturation pressure. Location of the steam saturation pressure eliminates ambiguity concerning the underlying phase state. For this fluid the spectral decomposition of the BIP matrix yields four nonzero eigenvalues. Since all eigenvalues are nonzero, the reduced variables formulation does not reduce the dimensionality of the linear system in the phase-split problem. In fact, the size of the linear system for three-phase flash is larger than that for the three-phase split problem in Case 1, despite the fluid in Case 1 consisting of twice the number of components. However, this is an unusual example, as a BIP matrix consisting of all nonzero interaction parameters is uncommon. The crux of the reduced variables method is the construction of an orthogonal basis from the eigenvalue decomposition of the BIP matrix. Case 4 examines the implications of this orthogonality for the solution of the nonlinear equations via the Newton method. Case 4: Four Pseudocomponent Crude Plus Water. Case 4 consists of the five component hydrocarbon−water mixture in Luo and Barrufet.64 Fluid parameters are summarized in Table 11. This is a pseudocomponent heavy crude. Binary interaction parameters provided by Zhu and

Figure 8. PT-phase diagram for the fluid in Case 4, taken from Luo and Barrufet.64 The mixture is comprised of 50% water and 50% hydrocarbon. Pressure is sampled every 0.25 bar, and temperature is sampled every 0.5 K.

We use this case study to emphasize the suitability of our reduced formulation for the solution of phase-split problems in which the presence of trace components adversely impacts the performance of Newton methods when using conventional variables. We perform phase-split calculations across the threephase region using the logarithm of K-values as primary variables (see Haugen et al.49 for details). We compute the condition number of the Jacobian matrix at convergence for both the conventional log K variables and the set of reduced variables. The condition number κ is computed as the ratio of the largest to smallest singular value, resulting from the singular value decomposition of a matrix: κ(J) =

mole fraction

MW (g/mol)

Tc (K)

Pc (bar)

ω

H2O PC1 PC2 PC3 PC4

0.50 0.15 0.10 0.10 0.15

18.015 30.00 156.00 310.00 400.00

647.3 305.556 638.889 788.889 838.889

220.8900 48.82 19.65 10.20 7.72

0.344 0.098 0.535 0.891 1.085

(36)

where J is the matrix and σ refers to a singular value. In our implementation we compute the condition number using the Eigen library.65 A comparison of Figures 9 and 10 shows that the reduced method yields a Jacobian that is far better conditioned. Both our reduced variables method implementation and our implementation of the log K method involve the simultaneous solution of all primary unknowns. However, the residual equations in the log K formulation exhibit extreme sensitivity to phase amounts. Decoupling of the Rachford−Rice equations from the fugacity residual

Table 11. Fluid Properties for Case 4, Taken from Luo and Barrufet64 component

σmax(J) σmin(J)

K


Article


Case 5: Nine-Component Condensate Plus Water. Case 5 combines water with the nine-component condensate from the SPE third comparative solution project.66 Fluid properties are provided in Table 13. Binary interaction parameters with water are estimated (see Table 14). For this fluid the spectral decomposition of the BIP matrix yields eight nonzero eigenvalues. Here we focus on convergence to the solution of the phase split problem. We use a mixture of 10% water and 90% hydrocarbon for our test case. A large three-phase region is present in the PTphase diagram in Figure 11. Figure 12 shows numerical convergence to the equilibrium solution at P = 200 bar, T = 290.5 K. In this example, successive substitution takes 2443 iterations to converge to the equilibrium result. On the other hand, the combined SS−Newton algorithm requires only five Newton iterations to converge. The switch criteria is set at 0.02, and 87 successive substitution iterations are required prior to switching to Newton. Results of the three-phase split are summarized in Table 15.

Figure 9. Condition number computed at the solution in the threephase region for the PT-phase diagram in Figure 8. Here, the condition number of the Jacobian matrix is associated with the log K formulation.

■

DISCUSSION Strongly nonideal behavior is characteristic of hydrocarbon− water systems and is particularly apparent in mixtures comprising a liquid hydrocarbon phase and an aqueous phase.51 Immiscibility is a presupposition in an isothermal context. Michelsen51 noted that phase-split calculations for these systems converge relatively quickly. In an isothermal setting a priori knowledge of water−oil immiscibility motivates preallocation of the phase state in lieu of performing stability analysis. Furthermore, phase-split calculations are particularly amenable to successive substitution in this context, as the Kvalues span several orders of magnitude. However, in the simulation of thermal recovery processes a wide range of operating conditions may be encountered, yielding several possible phase states. In addition, the high temperatures give rise to mutual solubility of water and hydrocarbons. In this research we have developed a strategy for reliable stability analysis and rapid phase-split calculations across the pressure− temperature parameter space. Although there is a large body of research dedicated to the resolution of phase behavior in hydrocarbon−water systems, the K-value approximation of phase behavior has persisted in the reservoir simulation community for several decades.67 Recently, however, there has been an increased interest in more accurately quantifying phase behavior in thermal recovery processes. One development is the use of free-water flash approximations of phase behavior.61,68,69 Another approach involves parametrization of compositional space.70

Figure 10. Condition number computed at the solution in the threephase region for the PT-phase diagram in Figure 8. The condition number corresponds to the Jacobian matrix associated with the reduced-variables formulation. The reduced Jacobian is far better conditioned.

equations eliminates this problem, although this is no longer a fully coupled solution, and the number of iterations required for convergence increases.22 It is possible that an alternative reduced variables method may yield comparable or even superior conditioning. However, the spectral decomposition is a very stable numerical algorithm due to the fact that it produces a unitary basis. Table 13. Component Properties for the Fluid in Case 566 component

mole fraction

MW (g/mol)

Tc (K)

Pc (bar)

ω

H2O CO2 N2 C1 C2 C3 C4−C6 PC1 PC2 PC3

0.10000 0.01089 0.01746 0.0.59391 0.07821 0.05319 0.08703 0.042705 0.013635 0.00297

18.015 44.01 28.013 16.043 30.07 44.097 66.86942 107.77943 198.56203 335.1979

647.3 304.7 126.2 190.6 305.43 369.8 448.08 465.62 587.8 717.72

220.8900 73.86796 33.94563 46.04085 48.83673 42.65743 35.50565 28.32348 17.06905 11.06196

0.344 0.225 0.04 0.013 0.0986 0.1524 0.21575 0.3123 0.5567 0.91692

L


Article

Industrial & Engineering Chemistry Research Table 14. Binary Interaction Parameters for the Fluid in Case 5a

a

BIP

H2O

CO2

N2

C1

C2

C3

C4−C6

PC1

PC2

PC3


0 0.0952 −0.48 −0.48 0.45 0.53 0.50 0.50 0.50 0.50

0 0 0 0.13 0.135 0.1277 0.1 0.1 0.1

0 0 0.05 0.08 0.1002 0.1 0.1 0.1

0 0 0 0.09281 0.130663 0.130663 0.130663

0 0 0 0.006 0.006 0.006

0 0 0.006 0.006 0.006

0 0 0 0

0 0 0

0 0

0

66

Hydrocarbon BIP data taken from Kenyon and Behie.

Hydrocarbon−water interaction parameters are estimated.

Table 15. Three-Phase Flash Results for the Fluid in Case 5a component

phase 1

phase 2

phase 3


0.000271381 0.0121475 0.0181017 0.630346 0.0889628 0.0624378 0.108759 0.0547045 0.0193717 0.00489795

0.000273447 0.0120496 0.0205266 0.685507 0.0850454 0.0561403 0.0860664 0.0410583 0.0114367 0.00189572

0.999972 1.24411e−05 1.48698e−05 3.65133e−07 2.49001e−10 2.59097e−14 1.11805e−18 3.89899e−24 8.29673e−49 2.76806e−90

a

Results obtained at P = 200 bar, T = 290.5 K. The hydrocarbon phases (phase 1 and phase 2) are near critical.

Figure 11. PT-phase diagram for mixture consisting of 10% water and 90% SPE3 condensate (Case 5). A large three-phase region is present. Pressure is sampled in increments of 0.25 bar, and temperature is sampled in increments of 0.25 K.

component partitioning across phases is not permitted. This is a good approximation in many cases, particularly at lower temperatures, although it is not completely general purpose. In the case of compositional space parametrization (CSP), the integrity of the parametrization is contingent on robust underlying phase equilibrium calculations with an equation of state. In addition, multiphase CSP assumes that two-phase tie-lines lie on the plane associated with a three-phase tiesimplex. However, tilted tie-lines may not align with the tiesimplex plane.71 Given this context, there is a need for thermodynamically rigorous quantification of hydrocarbon− water phase behavior in thermal reservoir simulation. Michelsen51 noted that difficult hydrocarbon−water systems are those which exhibit fairly narrow regions of immiscibility. Vigilant resolution of phase boundaries is rarely discussed in the literature, despite the potential implications of incorrect phase identification in reservoir simulation. Consistent navigation of phase boundaries is particularly important in thermal reservoir simulation in which the energy equation is tightly coupled to the phase state. Seemingly insignificant regions of immiscibility may be important in this context. This can be seen in comparison of Figure 4 and Figure 5. Failure to resolve the three phase state in pressure−temperature space is likely to result in a large error in the calculation of thermal properties, such as enthalpy or internal energy. Stability analysis carries the advantage of providing an initial estimate for compositions in phase-split calculations. Michelsen51 noted that this estimate is particularly accurate close to the phase boundary, where convergence of the phase-split calculation is generally most challenging. One of the key difficulties in stability analysis is the number of initial estimates

Figure 12. Analysis of convergence for a three-phase-split calculation at P = 200 bar and T = 290.5 K for the fluid in Case 5. The convergence criterion for the fugacity residual norm is set to 1e−10. Successive substitution requires 2443 iterations for convergence. This is attributable to the near criticality of the hydrocarbon phases. If we set the SS−Newton switch threshold to 0.02, we require 87 successive substitution iterations to reach the switch residual, and only 5 Newton iterations for convergence.

Both free-water flash and compositional space parametrization are improvements upon the legacy K-value approximation. Nevertheless, these recent advancements still constitute approximations of the underlying phase behavior. Free-water flash predetermines component insolubility, such that arbitrary M


Article


lations.40−43 It is not the purpose of this paper to provide a comprehensive comparison of conventional and reduced formulations for phase-split calculations. Such an analysis is difficult for two key reasons, as noted by Michelsen et al.41 First, an objective comparison of performance is made complicated by implementation differences. Second, reducedvariables formulations may result in the solution of a system that differs from the underlying full-dimensional problem. Certain reduced-variables formulations may create such differences via modification of the underlying mixing rules.33,34,36 Our formulation uses a linear combination of eigenvalues and eigenvectors from the BIP matrix of the energetic term. We retain all significant eigenvalues to ensure we faithfully represent the underlying physical system, with a cutoff absolute value of 1 × 10−14. In addition, we assess convergence by computing the residual of the fugacity norm, rather than the norm of a residual in the reduced parameter space.42 Finally, our formulation includes no modification of the mixing rules used in the PR EOS. We use the van Der Waals mixing rules with a single set of binary interaction parameters, as is standard practice in the reservoir simulation community. Hydrocarbon−water mixtures are characterized by the presence of trace components, which lead to ill-conditioned linear systems when Newton methods are used to solve the phase-split problem. The conventional approach to solving phase-split calculations consists of using independent variables such as component mole numbers, the logarithm of phase equilibrium ratios, or simply the equilibrium ratios themselves.76 In the case of component mole numbers, fugacity derivatives involving trace components lead to rank deficiency of the Jacobian or Hessian matrix. The presence of trace components in the aqueous phase leads to a near singular matrix. This was observed by Patterson et al.,77 who noted that the Jacobian became almost singular owing to compositions as low as 10−200. This observation was made for isenthalpic flash, but the rationale is equally applicable to the isothermal case. A solution to the numerical challenge of trace components is to use individually tailored variables for each component.51,78 However, this is cumbersome, and in the case of hydrocarbon−water mixtures at elevated temperatures there is no way of knowing how the components will partition a priori across indeterminate phases. In thermal reservoir simulation, high temperatures and pressures lead to mutual solubility of hydrocarbons and water. Another approach is empirical modification of the matrix, which may adversely affect convergence.78 Using the logarithm of K-values as primary variables avoids the problem of rank deficiency but results in a poorly conditioned linear system as derivatives of the fugacity eqs (eq 8 and eq 9) may be very large with respect to small changes in the logarithm of a phase-equilibrium ratio involving a trace component. The result of this poor conditioning is clearly shown in Figure 9. In contrast, our reduced-variables formulation results in a condition number that is 3 to 4 orders of magnitude smaller than that of the log K formulation in Figure 10. The reduced-variables formulation introduced in this paper is particularly well-suited to hydrocarbon−water systems. The advantage of constructing a parameter space using a spectral decomposition of the BIP matrix is twofold: (i) the resulting reduced Jacobian matrix is full rank, by-construction; (ii) the size of the linear system is substantially reduced, provided there is a significant number of zero values in the BIP matrix.

required. Standard practices of estimating K-values for hydrocarbon systems do not work for hydrocarbon−water mixtures. Our approach to partitioning the domain using the steam saturation curve allows us to identify how many initial K-value estimates are required a priori. The use of multiple two-phase flash calculations above the steam saturation pressure enables us to reliably identify the correct phase state and provide a good initial estimate for ensuing flash calculations. Our solution strategy addresses the difficulties encountered with narrow regions of three-phase behavior and multiple two-phase states. Case 1 shows that in these scenarios, the presence of a three-phase state may not be predicted by performing stability analysis using the minimum Gibbs free energy two-phase flash result. Initiating phase-split calculations in reservoir simulation without prior stability analysis may be particularly error prone for fluids characterized by capricious regions of three phase behavior, as in Cases 1, 2, and 3. We are the first to propose a systematic stability analysis strategy for a wide range of temperatures and pressures likely to be encountered in simulation of thermal recovery. Several authors have noted that test phase selection does not influence the outcome of stability testing in multiphase mixtures.22,72,73 The reasoning behind this assertion is that the outcome of stability testing using different trial phases will be the same, based on Gibbs free energy analysis of the equilibrium phases.45 However, the tangent plane criterion only guarantees phase instability in the event a negative TPD is detected. Stability test calculations are performed using locally convergent methods. Given the same initial 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. The decision to perform stability testing with one phase instead of another may render the stability analysis far more difficult and far less likely to detect a stationary point with a negative TPD. For this reason alone, it is imperative that stability testing is guided by the identity and composition of the phases present following two-phase flash. For hydrocarbon−water mixtures, phases are readily identifiable. In compositional reservoir simulation, there exists a nonlinear coupling of the local thermodynamic constraints to the global solution of conservation laws.74,75 In thermal compositional simulation, the energy equation exhibits extreme sensitivity to phase behavior at each grid block. For this reason, the accurate solution of phase-split problems is important to the fidelity of the global Newton loop. This is especially the case when molar variables are used in the solution of the conservation laws. In this case, convergence of local thermodynamic constraints at the grid block level is required to close the system of equations. This cannot always be achieved using successive substitution alone, and a secondorder method is required. Case 5 is a relevant example. In a three-phase-split calculation the vapor and liquid hydrocarbon phases may approach critical behavior. Figure 12 clearly shows that as bicritical behavior is approached for a hydrocarbon− water mixture, successive substitution does not converge to the required tolerance in a reasonable number of iterations. The literature on phase equilibrium calculations contains a great deal of conjecture regarding the efficacy of reducedvariables formulations for the solution of phase-split calcuN


Article


Figure A1. Flow chart of the hydrocarbon−water mixture stability analysis.

In regards to the first point, the set of reduced-variables acts as a low-pass filter on the underlying physics of the component interactions manifest in the BIP matrix. The eigenvectors obtained from the spectral decomposition of the BIP matrix are linearly independent, even in the case of trace components.

As a result of this, the standard SSI−Newton solution protocol can be used with confidence for phase-split calculations. The linearized system can be reduced in size because most hydrocarbon−water mixtures consist of a large number of zero BIP terms. For this reason, the reduced-variables method O


Article

Industrial & Engineering Chemistry Research can be used to solve phase-split calculations for complex heavy oils with a a relatively small system. The size of the linear system scales with the rank of the BIP matrix. As a result, the reduced basis facilitates more detailed fluid characterization, in which components are lumped into a larger number of pseudocomponents. In turn, improved fluid characterization may enable more faithful representation of complex thermal recovery mechanisms in thermal EOR techniques such as SAGD and enhanced-SAGD.

∑

G̅ =

Nc

niV ln(yi )

+

i=1

i=1

Nc

+

∑ niW ln(wi)

(B-2)

i=1

G̅ E = Vg ̅ E,V + Lg ̅ E,L + Wg ̅ E,W E,V

E,L

(B-3)

E,W

where g̅ , g̅ , and g̅ are the excess Gibbs free energies of the vapor, liquid, and aqueous phases, respectively. Thus

■

Nc

G̅ =

Nc

Nc

∑ niV ln(yi ) + ∑ niL ln(xi) + ∑ niW ln(wi) + Vg ̅ E,V i=1

i=1

+ Lg ̅

E,L

+ Wg ̅

i=1

E,W

(B-4)

However, the total moles of the feed is the sum of the moles across the phases (B-5)

V+L+W=1

Thus, if QV, QW, V, and W are independent variables QL =

Q F − VQV − WQW L

(B-6)

where the superscript F refers to moles of the feed and n L = n − n V − nW

(B-7)

Thus, in functional form n = (n1 , n2 , ... ., n Nc)

(B-8)

nV = (n1V , n2V , ..., n NVc)

(B-9)

n L = (n1L , n2L , ..., n NLc)

(B-10)

nW = (n1W , n2W , ..., n NWc )

(B-11)

Equations B-5, B-6, and B-7 in eq B-4 yield Nc

G̅ =

Nc

Nc

∑ niV ln(yi ) + ∑ niW ln(wi) + ∑ (ni − niV i=1

i=1

− niW )ln(xi)

i=1

+ Vg ̅ E,V (Q V ) + Wg ̅ E,W (QW ) + (F − V − W )g ̅ E,L (Q V , QW )

(B-12)

Via the definition Q̅ V = V QV and Q̅ W = WQW we obtain

■

Nc

G̅ =

APPENDIX A: HYDROCARBON−WATER MIXTURE STABILITY ANALYSIS A flow chart for this analysis is shown in Figure A1.

Nc

Nc

∑ niV ln(yi ) + ∑ niW ln(wi) + ∑ (ni − niV

− niW )ln(xi)

Wy i V ij Q̅ V yz i Wy zz + Wg E,W jjj Q̅ zzz + (F − V − W )g E,L jjj Q̅ , Q̅ zzz + Vg ̅ E,V jjjj zz zz ̅ jj W zz ̅ jj V V W { k { k { k i=1

i=1

i=1

(B-13)

APPENDIX B: MATHEMATICAL DEVELOPMENT OF THE REDUCED METHOD FOR THREE-PHASE FLASH

The expression is similar to the form provided by Petitfrere and Nichita60 for multiphase Gibbs free energy with reduced variables. To be clear, the component mole numbers nVi and nW i are dependent variables. The variables are subject to the following constraints:

Development of the Reduced Variables

The Gibbs free energy of mixture is the summation of the ideal Gibbs free energy and the excess Gibbs free energy: G̅ = G̅ I + G̅ E

∑

niL ln(xi)

The excess Gibbs energy of the mixture is given by

CONCLUSIONS We have developed a strategy for phase-equilibrium calculations involving hydrocarbon−water mixtures. Our approach uses a cubic EOS without modification of the van Der Waals mixing rules. The strategy is sensitive to the distinctive properties of hydrocarbon−water mixtures, including the physical and numerical implications of a near-pure-component aqueous phase. For stability analysis, we use the steam saturation curve to partition the pressure temperature domain, and we use a maximum of three sets of initial K-value estimates to initiate each stability test. In addition, we perform multiple two-phase flash calculations to consistently resolve phase boundaries. We have verified that this approach is robust through rigorous testing of a variety of fluids from the literature. To the best of our knowledge, we are the first to use the steam saturation curve to guide the initial estimate of K-values in stability testing. We introduce a three-phase reduced-variables formulation to address shortcomings with existing phase-split calculations for hydrocarbon−water mixtures. We have demonstrated excellent performance of the reduced formulation for several hydrocarbon−water mixtures. The set of reduced-variables yields a full-rank Jacobian matrix when Newton methods are used to solve the phase-split problem. In comparison with the alternative formulation, the natural logarithm of K-values, the reduced method yields a Jacobian matrix with a condition number several orders of magnitude lower. We expect the developments in this paper to have implications for phase-equilibrium computations in a nonisothermal setting (isenthalpic, isentropic). In addition, we believe this research serves as an important step toward strengthening the coupling of phase-equilibrium calculations and the transport equations in thermal compositional simulation.

■

Nc

I

Q̅ αV =

(B-1)

Nc

∑ qα niV α = 1, ..., m i=1

Here, we work in terms of normalized Gibbs free energy, G̅ = G/RT. For a three-phase VLW system the ideal Gibbs free energy is the combined ideal energies of the three phases:

Q̅ αW =

(B-14)

Nc

∑ qα niW α = 1, ..., m i=1

P

i

i

(B-15) DOI: 10.1021/acs.iecr.9b00695 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Nc

V

∑ yb i i

b =

bW =

Nc

(B-17)

niV

∑

∑ wbi i − bW

Nc

(B-19)

i=1

−

i=1

Q αV

= Rα

Jacobian Matrix

Here we provide the full form of the Jacobian matrix for the three-phase split calculation using the Newton method with reduced variables. We explicitly include the Jacobian matrix to facilitate reproducibility. Okuno56 provides a detailed description of Jacobian elements for an alternate three-phase reduced variables method.

α = 1, ..., m

= RM

ÅÄÅ ∂R ÅÅ 1 ÅÅ ÅÅ ∂Q V ÅÅ 1 ÅÅ ÅÅ ÅÅ ∂ ÅÅ ÅÅ ∂R ÅÅ M + 1 ÅÅ ÅÅ ∂Q V ÅÅ 1 J = ÅÅÅ ÅÅ ∂R ÅÅ M + 2 ÅÅ ÅÅ ∂Q V 1 ÅÅ ÅÅ ÅÅ ∂ ÅÅ ÅÅ ÅÅ ∂R ÅÅ 2M + 2 ÅÅ ÅÅ ∂Q V ÅÇ 1

(B-21)

μ μ μ μ

∂R1

∂R1

∂Q mV

∂bV

μ ∂RM + 1

μ ∂RM + 1

∂Q mV

∂bV

∂RM + 2

∂RM + 2

∂Q mV

∂R1

∂R1 ∂V

V

∂b

∂Q γV ∂R k V

∂b

=

∑ q αi i=1 Nc

=

∂R k = ∂V

∑ q αi i=1 Nc

∑ q αi i=1

∂Q γV

∂R1

∂Q mW

∂bW

μ μ μ ∂RM + 1 ∂RM + 1 μ W ∂Q 1 ∂Q mW

μ ∂RM + 1

∂RM + 2 ∂V

∂RM + 2

∂RM + 2

∂RM + 2

∂Q mW

∂bW

μ

∂Q 1W

μ

∂bW

μ μ μ μ μ μ μ ∂R 2M + 2 ∂R 2M + 2 ∂R 2M + 2 ∂R 2M + 2 ∂R 2M + 2 ∂R 2M + 2 μ μ W V V ∂ V ∂Q 1 ∂Q mW ∂bW ∂Q m ∂b Nc

∂R k

=

∂Q γW ∂R k ∂bW

− δαγ (B-26)

∑ q αi i=1 Nc

Nc i=1

∂yi ∂Q γW

∂bW

i=1

∑ q αi

ÑÉÑ ÑÑ ÑÑ ÑÑ ÑÑ ÑÑ Ñ μ ÑÑÑÑ Ñ ∂RM + 1 ÑÑÑÑ ÑÑ ∂W ÑÑÑ ÑÑ ÑÑ ∂RM + 2 ÑÑÑ ÑÑ ∂W ÑÑÑÑ ÑÑ Ñ μ ÑÑÑÑ ÑÑ ∂R 2M + 2 ÑÑÑÑ Ñ ∂W ÑÑÑÑ Ö ∂R1 ∂W

(B-29)

∂yi

∑ q αi

=

∂R k = ∂W

(B-30)

∂yi ∂W

(B-31)

where k, α, γ = 1, 2, ..., m. For the derivatives of the residual RM

∂yi

∂RM

V

∂Q γV

∂b

∂R1

μ ∂RM + 1 ∂V

Here we give the expressions for the partial derivatives of residuals R1−RM+1 (eqs B-21−B-22) with respect to the reduced variables. The partial derivatives of residuals RM+2− R2M+2 (eqsB-23−B-25) are analogous. We begin with the derivatives of residuals R1−Rm. ∂yi

μ

∂Q 1W

Partial Derivatives Constituting the Jacobian Matrix eElements

Nc

(B-25)

In eqs B-21−B-25 we take M = m + 1. Note that R is used here to refer to the residuals of the nonlinear equations and is unrelated to the gas constant.

(B-20)

i=1

∂R k

(B-24)

i=1

Nc

− bV ∑ yb i i

= R 2M + 1

Nc

The phase-split calculations are performed via solution of the system of nonlinear equations. Here we outline the procedure via (i) successive substitution and (ii) Newton’s method. The independent set of reduced variables is {QV, bV, V, QV, bW, W}. The nonlinear equations for three-phase flash are as follows: i

(B-23)

∑ (wi − xi) = R 2M + 2

Solution of the Nonlinear Equations

Nc

α = 1, ..., m

i=1

∑ niW

∑ yq i α

i

Nc

(B-18)

i=1

− Q αW = R α + M + 1

∑ wqi α i=1

Nc

W=

(B-22)

i=1

Nc

∑ wbi i i=1

V=

∑ (yi − xi) = RM + 1

(B-16)

i=1

(B-27)

∂yi

∂RM

∂V

∂Q γW

(B-28) Q

Nc

=

∑ bi i=1

∂yi ∂Q γV

Nc

=

∑ bi i=1

(B-32)

∂yi ∂Q γW

(B-33) DOI: 10.1021/acs.iecr.9b00695 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


∂RM

=

V

∂b

∑ bi i=1 Nc

∂RM

=

∂bW

∑ bi i=1

∂RM = ∂V ∂RM = ∂W

Nc

∑ bi i=1 Nc

∑ bi i=1

∂yi

Of course, yi can be written as

−1

V

∂b

(B-34)

∂bW

∂yi

∂yi ∂V

∂yi

∂yi

∂Q αW

(B-37)

ij ∂y ∂xi yzzz j zz = ∑ jjj i V − jj ∂Q ∂Q γV zz i=1 k γ {

∂yi ∂b

∂Q γW ∂RM + 1

=

∂yi

(B-38)

∂b

∂b

i=1

=

W

∂b

i=1

∂RM + 1 = ∂V ∂RM + 1 = ∂W

k ∂b k ∂b

−

V

i ∂yi

W

i ∂yi

∑ jjjjj Nc

i=1

k ∂V

k ∂W

−

i ∂yi

i=1

(B-39)

∂V

(B-40)

∂xi yzz zz ∂bW z{

∂W

∂xi V

∂b

∂xi ∂bW ∂xi ∂V ∂xi ∂W

∂Q αW

V

∂b

∂K iV W

∂b

∂Q αV

+ K iV

+ K iV + K iV

(B-51)

∂xi ∂Q αW

(B-52)

∂xi ∂bV

(B-53)

∂xi ∂bW

(B-54)

(B-41)

(B-42)

∂xi yzz z ∂W zz{

∂K iV

(B-43)

i ∂K V ∂K W y −zijjjV iW + W iW zzz ∂ ∂ Q Qα { α k = V [1 + V (K i − 1) + W (K iW − 1)]2

(

− zi V =

[1 +

V (K iV

(

∂bV

+W

− 1) +

− zi V =

∂K iV

∂K iV W

∂b

∂K iW ∂bV

∂K iV

+W

∂K iW ∂bW

∂bV

∂K iV

(B-44)

∂bW

(B-45)

)

W (K iW

− 1)]2

(B-46)

)

[1 + V (K iV − 1) + W (K iW − 1)]2 ÄÅ É Å Ñ ∂K V ∂K W Ñ ziÅÅÅÅ1 − K iV − V ∂Vi − W ∂Vi ÑÑÑÑ ÅÇ ÑÖ = [1 + V (K iV − 1) + W (K iW − 1)]2 É ÅÄÅ Ñ ∂K V ∂K W Ñ ziÅÅÅÅ1 − K iW − V ∂Wi − W ∂Wi ÑÑÑÑ Å ÑÖ Ç = V W [1 + V (K i − 1) + W (K i − 1)]2

(B-55)

∂K iV ∂x + K iV i ∂W ∂W

(B-56)

= xi

∂Q αW

i ∂K V ∂K W y −zijjjV iV + W i V zzz ∂ ∂ Q Qα { α k = V [1 + V (K i − 1) + W (K iW − 1)]2

∂K iV ∂x + K iV i ∂V ∂V

Now, for the partial derivatives involving equilibrium ratios ÄÅ É VÑ L ÅÅ ÑÑ ̂ ̂ ∂ ϕ ∂ ϕ ln ln Å ∂K iV V ÑÑ i Ñ i VÅ Å = − + K Å i Å ÑÑ V Ñ ÅÅ L ∂Q L ∂Q αV ∂ Q ÅÅÇ ÑÑÖ α α Ñ (B-57)

Construction of Jacobian Elements

∂Q αW

= xi

∂yi

∂xi zyz z ∂V zz{

−

∑ jjjjj Nc

∂yi

∂xi zyz zz ∂bV z{

−

∑ jjjjj Nc

∂RM + 1

i ∂yi

∑ jjjjj Nc

V

∂xi

∂K iV

∂xi

Nc

∂RM + 1

∂Q αV

+ K iV

∂Q αV

∂K iV

= xi

W

∂K iV

= xi

= xi

V

ij ∂y ∂xi yzzz j = ∑ jjj iW − zz jj ∂Q ∂Q γW zz γ i=1 k {

∂xi

= xi

∂Q αV

(B-36)

∂W

(50)

Thus, we can write

(B-35)

Nc

∂Q γV

1 + V (K iV − 1) + W (K iW − 1)

∂yi

Now, for the derivatives of the residual RM+1 ∂RM + 1

K iV zi

yi = K iV xi =

■

(B-47)

=

W −K iL

∂ln ϕî

L

L ∂Q αL

ÄÅ ÉÑ ÅÅ ̂L ̂V ÑÑÑ ∂ ϕ ∂ ϕ ln ln Å V Å i Ñ i Ñ = −K iV ÅÅÅ + ÑÑ V Ñ ÅÅ L ∂bL ∂ b ÑÑ ÅÅÇ ÑÖ

= −K iL

W ∂ln ϕî L ∂bL

(B-58)

(B-59)

L

É ÅÄÅ m L L Ñ Ñ ∂ln(ϕî ) ∂ln(ϕî ) ÑÑÑ ∂K iV K iV ÅÅÅÅ L V L V ÑÑ ÅÅ ∑ (Q α − Q α ) = ( b b ) + − Ñ ∂V L ÅÅÅ α = 1 ∂Q αL ∂b L ÑÑÑ ÅÇ ÑÖ (B-61) É ÅÄÅ m L L Ñ Ñ ∂ln(ϕî ) ∂ln(ϕî ) ÑÑÑ ∂K iV K iV ÅÅÅÅ L W L W ÑÑ ÅÅ ∑ (Q α − Q α ) = ( b b ) + − Ñ ∂W L ÅÅÅ α = 1 ∂Q αL ∂b L ÑÑÑ ÅÇ ÑÖ (B-62) (B-60)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

(B-48)

Michael Connolly: 0000-0003-3962-9001 Huanquan Pan: 0000-0003-4229-971X Notes

The authors declare no competing financial interest.

(B-49) R


Article


■

(8) Zaydullin, R.; Voskov, D.; Tchelepi, H. Comparison of EoSBased and K-Values-Based Methods for Three-Phase Thermal Simulation. Transp. Porous Media 2017, 116, 663−686. (9) Johns, R.; Dindoruk, B.; Orr, F., Jr Analytical theory of combined condensing/vaporizing gas drives. SPE Advanced Technology Series 1993, 1, 7−16. (10) Pan, H.; Connolly, M.; Tchelepi, H. Multiphase equilibrium calculation framework for compositional simulation of CO2 injection in low temperature reservoirs. Ind. Eng. Chem. Res. 2019, 58, 2052− 2070. (11) Peng, D.-Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (12) Robinson, D. B.; Peng, D.-Y. The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs; Gas Processors Association, 1978. (13) Soave, G. Equilibrium constants from a modified RedlichKwong equation of state. Chem. Eng. Sci. 1972, 27, 1197−1203. (14) Economou, I. G.; Tsonopoulos, C. Associating models and mixing rules in equations of state for water/hydrocarbon mixtures. Chem. Eng. Sci. 1997, 52, 511−525. (15) Wu, J.; Prausnitz, J. M. Phase equilibria for systems containing hydrocarbons, water, and salt: An extended Peng- Robinson equation of state. Ind. Eng. Chem. Res. 1998, 37, 1634−1643. (16) Moortgat, J.; Li, Z.; Firoozabadi, A. Three-phase compositional modeling of CO2 injection by higher-order finite element methods with CPA equation of state for aqueous phase. Water Resour. Res. 2012, 48. DOI: 10.1029/2011WR011736. (17) Michelsen, M. L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1−19. (18) Sun, A. C.; Seider, W. D. Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy. Fluid Phase Equilib. 1995, 103, 213−249. (19) McDonald, C. M.; Floudas, C. A. Global optimization for the phase stability problem. AIChE J. 1995, 41, 1798−1814. (20) 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. (21) Harding, S.; Floudas, C. Phase stability with cubic equations of state: Global optimization approach. AIChE J. 2000, 46, 1422−1440. (22) Li, Z.; Firoozabadi, A. General strategy for stability testing and phase-split calculation in two and three phases. SPE Journal 2012, 17, 1−096. (23) Paterson, D.; Yan, W.; Michelsen, M.; Stenby, E. Isenthalpic Flash for Thermal Recovery; ECMOR XV-15th European Conference on the Mathematics of Oil Recovery, 2016. (24) Mahabadian, M. A.; Chapoy, A.; Burgass, R.; Tohidi, B. Development of a multiphase flash in presence of hydrates: Experimental measurements and validation with the CPA equation of state. Fluid Phase Equilib. 2016, 414, 117−132. (25) Pang, W.; Li, H. A. Application of augmented free-water Rachford-Rice algorithm to water/hydrocarbons mixtures considering the dissolution of methane in the aqueous phase. Fluid Phase Equilib. 2018, 460, 75−84. (26) Sabet, N.; Gahrooei, H. R. E. A new robust stability algorithm for three phase flash calculations in presence of water. J. Nat. Gas Sci. Eng. 2016, 35, 382−391. (27) Michelsen, M. L.; Heidemann, R. A. Calculation of critical points from cubic two-constant equations of state. AIChE J. 1981, 27, 521−523. (28) Michelsen, M. L. Simplified flash calculations for cubic equations of state. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 184− 188. (29) Hendriks, E. M. Reduction theorem for phase equilibrium problems. Ind. Eng. Chem. Res. 1988, 27, 1728−1732. (30) Hendriks, E.; Van Bergen, A. Application of a reduction method to phase equilibria calculations. Fluid Phase Equilib. 1992, 74, 17−34.

ACKNOWLEDGMENTS The authors thank Total S.A. for supporting this work through the Stanford Total Enhanced Modeling of Source Rock (STEMS) project. In addition, the authors acknowledge the support of the Stanford University Petroleum Research Institute for reservoir simulation (SUPRI-B).

■

G̅ G̅ I g̅E,j β ϕ̂ ji μ ω a b fjî G H KVi KW i L Nc nji Np P Pc Qk R T Tc V v W wi X xi yi Z zi

■

NOMENCLATURE normalized Gibbs free energy ideal normalized Gibbs free energy excess normalized Gibbs free energy of phase j matrix of binary interaction parameters fugacity coefficient for component i in phase j chemical potential acentric factor PR EOS energetic term PR EOS covolume term fugacity of component i in phase j Gibbs free energy enthalpy, J equilibrium constant for oil−vapor equilibrium constant for oil−water oleic or hydrocarbon liquid phase amount number of components component mole number for component i in phase j number of phases pressure, bar critical pressure, bar reduced variable k resulting from spectral decomposition of BIP matrix gas constant, 8.314 J/mol·K temperature, K critical temperature, K vapor phase amount volume, L aqueous phase amount mole fraction of component i in aqueous phase set of reduced variables mole fraction of component i in oleic phase mole fraction of component i in vapor phase compressibility factor mole fraction of component i in overall mixture REFERENCES

(1) Wilson, G. M. A modified Redlich-Kwong equation of state, application to general physical data calculations; 65th National AIChE Meeting, Cleveland, OH, 1969. (2) Ghosh, S.; Johns, R. T. An equation-of-state model to predict surfactant/oil/brine-phase behavior. SPE Journal 2016, 21, 1106− 1125. (3) Zhu, D.; Okuno, R. Robust Isenthalpic Flash for Multiphase Water/Hydrocarbon Mixtures. SPE Journal 2015, 20, 1350−1365. (4) Sun, H.; Raghavan, A.; Haugen, K.; Chempath, S. An improved isenthalpic flash algorithm based on maximization of entropy. Fluid Phase Equilib. 2017, 438, 18−28. (5) Orr, F. M. Theory of gas injection processes; Tie-Line Publications, 2007. (6) Varavei, A.; Sepehrnoori, K. An EOS-based compositional thermal reservoir simulator; SPE Reservoir Simulation Symposium, 2-4 February, 2009, The Woodlands, Texas; Society of Petroleum Engineers, 2009. (7) Heidari, M. Equation of State Based Thermal Compositional Reservoir Simulator for Hybrid Solvent/Thermal Processes. Ph.D. Thesis, University of Calgary, 2014. S


Article


(56) Okuno, R. Modeling of multiphase behavior for gas flooding simulation. Ph.D. Thesis, The University of Texas at Austin, 2009. (57) Okuno, R.; Johns, R.; Sepehrnoori, K. A new algorithm for Rachford-Rice for multiphase compositional simulation. SPE Journal 2010, 15, 313−325. (58) Younis, R. M. Modern advances in software and solution algorithms for reservoir simulation. Ph.D. Thesis, The University of Texas at Austin, 2011. (59) Zhou, Y. Parallel general-purpose reservoir simulation with coupled reservoir models and multisegment wells. Ph.D. Thesis, Stanford University, 2012. (60) Petitfrere, M.; Nichita, D. V. Multiphase equilibrium calculations using a reduction method. Fluid Phase Equilib. 2015, 401, 110−126. (61) Lapene, A.; Nichita, D. V.; Debenest, G.; Quintard, M. Threephase free-water flash calculations using a new Modified Rachford− Rice equation. Fluid Phase Equilib. 2010, 297, 121−128. (62) Brantferger, K. M. Development of a thermodynamically consistent, fully implicit, equation-of-state, compositional steamflood simulator. Ph.D. Thesis, The University of Texas at Austin, 1991. (63) Varavei, A. Development of an equation-of-state thermal flooding simulator. Ph.D. Thesis, The University of Texas at Austin, 2009. (64) Luo, S.; Barrufet, M. A. Reservoir simulation study of water-inoil solubility effect on oil recovery in steam injection process. SPE Reservoir Evaluation & Engineering 2005, 8, 528−533. (65) Guennebaud, G.; Jacob, B. Eigen v3. http://eigen.tuxfamily.org, 2010. (66) Kenyon, D.; Behie, G. Third SPE Comparative Solution Project: Gas Cycling of Retrograde Condensate Reservoirs. JPT, J. Pet. Technol. 1987, 39 (8), 981−997. (67) Aziz, K.; Ramesh, A.; Woo, P. Fourth SPE comparative solution project: comparison of steam injection simulators. JPT, J. Pet. Technol. 1987, 39, 1576−1584. (68) Tang, Y.; Saha, S. An Efficient Method to Calculate ThreePhase Free-Water Flash for Water−Hydrocarbon Systems. Ind. Eng. Chem. Res. 2003, 42, 189−197. (69) Heidari, M.; Nghiem, L. X.; Maini, B. B. Improved Isenthalpic Multiphase Flash Calculations for Thermal Compositional Simulators; SPE Heavy Oil Conference, Canada; Society of Petroleum Engineers, 2014. (70) Iranshahr, A.; Voskov, D.; Tchelepi, H. Tie-simplex parameterization for EoS-based thermal compositional simulation. SPE Journal 2010, 15, 545−556. (71) Zaydullin, R. Compositional space parameterization methods for thermal-compositional simulation. Ph.D. Thesis, Stanford University, 2015. (72) Nghiem, L. X.; Li, Y.-K. Computation of multiphase equilibrium phenomena with an equation of state. Fluid Phase Equilib. 1984, 17, 77−95. (73) Perschke, D. R. Equation-of-state phase-behavior modeling for compositional simulation. Ph.D. Thesis, The University of Texas at Austin, 1988. (74) Young, L. C.; Stephenson, R. E. A generalized compositional approach for reservoir simulation. SPEJ, Soc. Pet. Eng. J. 1983, 23, 727−742. (75) Acs, G.; Doleschall, S.; Farkas, E. General purpose compositional model. SPEJ, Soc. Pet. Eng. J. 1985, 25, 543−553. (76) Nghiem, L.; Heidemann, R. General acceleration procedure for multiphase flash calculation with application to oil-gas-water systems; Proceedings of the 2nd European Symposium on Enhanced Oil Recovery, 1982; pp 303−316. (77) Paterson, D.; Yan, W.; Michelsen, M. L.; Stenby, E. H. Robust and Efficient Isenthalpic Flash Algorithms for Thermal Recovery of Heavy Oil; SPE Improved Oil Recovery Conference; Society of Petroleum Engineers, 2016. (78) Michelsen, M. L.; Mollerup, J.; Breil, M. P. Thermodynamic Models: Fundamental & Computational Aspects; Tie-Line Publications, 2008.

(31) Firoozabadi, A.; Pan, H. Fast and Robust Algorithm for Compositional Modeling: Part I-Stability Analysis Testing. SPE Journal 2002, 7, 78−89. (32) Pan, H.; Firoozabadi, A. Fast and Robust Algorithm for Compositional Modeling: Part II-Two-Phase Flash Computations. SPE Journal 2003, 8, 380−391. (33) Li, Y.; Johns, R. T. Rapid flash calculations for compositional simulation. SPE Reservoir Evaluation & Engineering 2006, 9, 521−529. (34) Okuno, R.; Johns, R. T.; Sepehrnoori, K. Three-phase flash in compositional simulation using a reduced method. SPE Journal 2010, 15, 689−703. (35) Nichita, D. V.; Graciaa, A. A new reduction method for phase equilibrium calculations. Fluid Phase Equilib. 2011, 302, 226−233. (36) Gorucu, S. E.; Johns, R. T. New reduced parameters for flash calculations based on two-parameter BIP formula. J. Pet. Sci. Eng. 2014, 116, 50−58. (37) Nichita, D. V.; Broseta, D.; de Hemptinne, J.-C. Multiphase equilibrium calculation using reduced variables. Fluid Phase Equilib. 2006, 246, 15−27. (38) Saaf, F. Efficient Application Of Reduced Variable Transformation And Conditional Stability Testing In Reservoir Simulation F l as h C a l c u l at i o n s . h t t p s : / / w w w . g o o g l e . co m / p at e nt s / US20070282582, US Patent App. US11/758,215, 2007. (39) Schlumberger. INTERSECT: Technical Description 2015.1, 2015. (40) Haugen, K. B.; Beckner, B. A critical comparison of reduced and conventional eos algorithms. SPE Journal 2013, 18, 378−388. (41) Michelsen, M.; Yan, W.; Stenby, E. H. A Comparative Study of Reduced-Variables-Based Flash and Conventional Flash. SPE Journal 2013, 18, 952−959. (42) Gorucu, S. E.; Johns, R. T. Comparison of reduced and conventional two-phase flash calculations. SPE Journal 2015, 20, 294−305. (43) Petitfrere, M.; Nichita, D. V. A comparison of conventional and reduction approaches for phase equilibrium calculations. Fluid Phase Equilib. 2015, 386, 30−46. (44) Mohebbinia, S.; Sepehrnoori, K.; Johns, R. T. Four-phase equilibrium calculations of carbon dioxide/hydrocarbon/water systems with a reduced method. SPE Journal 2013, 18, 943−951. (45) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs energy analysis of phase equilibria. SPEJ, Soc. Pet. Eng. J. 1982, 22, 731−742. (46) Rasmussen, C. P.; Krejbjerg, K.; Michelsen, M. L.; Bjurstrøm, K. E. Increasing the Computational Speed of Flash Calculations With Applications for Compositional, Transient Simulations. SPE Reservoir Evaluation & Engineering 2006, 9, 32−38. (47) Pan, H.; Tchelepi, H. A. Compositional flow simulation using reduced-variables and stability-analysis bypassing; SPE Reservoir Simulation Symposium; Society of Petroleum Engineers, 2011. (48) Zaydullin, R.; Voskov, D.; Tchelepi, H. Phase-state identification bypass method for three-phase thermal compositional simulation. Computational Geosciences 2016, 20, 461−474. (49) Haugen, K. B.; Firoozabadi, A.; Sun, L. Efficient and robust three-phase split computations. AIChE J. 2011, 57, 2555−2565. (50) Rachford, H., Jr; Rice, J. Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium. JPT, J. Pet. Technol. 1952, 4, 19−3. (51) Michelsen, M. L. The isothermal flash problem. Part II. Phasesplit calculation. Fluid Phase Equilib. 1982, 9, 21−40. (52) Gorucu, S.; Johns, R. T. Robustness of three-phase equilibrium calculations. J. Pet. Sci. Eng. 2016, 143, 72−85. (53) Saul, A.; Wagner, W. International equations for the saturation properties of ordinary water substance. J. Phys. Chem. Ref. Data 1987, 16, 893−901. (54) Eubank, P. T.; Wu, C. H.; Alvarado, J. F.; Forero, A.; Beladi, M. K. Measurement and prediction of three-phase water/hydrocarbon equilibria. Fluid Phase Equilib. 1994, 102, 181−203. (55) Michelsen, M. Phase equilibrium calculations. What is easy and what is difficult? Comput. Chem. Eng. 1993, 17, 431−439. T


Article

Industrial & Engineering Chemistry Research (79) Venkatramani, A. V.; Okuno, R. Modeling of Multiphase Behavior for Water/n-Alkane Mixtures by Use of the Peng-Robinson EOS; SPE Heavy Oil ConferenceCanada; Society of Petroleum Engineers, 2014.

U


Three-Phase Equilibrium Computations for HydrocarbonâWater

Recommend Documents

Three-Phase Equilibrium Computations for HydrocarbonâWater