An improved reduction method for phase stability testing in the single-phase region

1 Introduction

Phase stability analysis is very important in process systems engineering and petroleum and gas reservoir, production and transport engineering. Phase stability testing assesses the state of a mixture at given specifications and is essential in the initialization of flash calculations, in checking the results of phase split calculations and in phase diagram construction. It consists in an unconstrained minimization of the tangent plane distance (TPD) function [1]. The stability test is performed in most cases at a given pressure and temperature, but many implementations using various other specifications were proposed [2,3].

The stability test limit locus (STLL) [45,6,7,] is an important underlying property of multicomponent system phase diagrams, because in its vicinity, the number of iterations for phase stability testing dramatically increases and divergence may occur. The Hessian matrix is singular at the STLL for a nontrivial trial phase composition. The cause of convergence problems in phase stability calculations in the single-phase state is the topology of the TPD surface, as explained in refs. [6,7]. At the STLL, the TPD surface exhibits a saddle point, and for conditions outside the STLL in the pressure–temperature plane at a given composition, the iterates have to cross a domain of indefiniteness of the Hessian matrix, starting from one of the two-sided initial guesses [7]. This makes stability testing in the vicinity of the STLL really challenging, and any gradient-based algorithm will experience difficulties in this region [56,7,8,].

The reduction methods [9,10], in which the solution is sought in the hyperspace defined by the reduction parameters, rather than in the compositional hyperspace, are very attractive alternatives to the conventional methods, especially for mixtures with many components, which can often be encountered in refining and chemical process simulators (mixtures may have hundreds of components, resulting from mixing several feeds) or in petroleum engineering applications. In the reduction methods, the number of independent variables does not depend on the number of components in the mixture but only on the number of components having nonzero binary interaction parameters (BIPs) in the equation of state (EoS) with the remaining ones. Reduced stability testing routines were included in several commercial and academic compositional reservoir simulators, in which thermodynamic equilibrium equations are coupled with equations of fluids flow through porous media and need to be solved a huge number of times for full-field simulations.

In this work, a new reduction method (which is an improved version of our previous formulation) for the phase stability problem is presented, consisting in Newton iterations (with a switch-back to successive substitution iterations [SSI], as a safety feature) with a new set of independent variables and residual functions. The paper is structured as follows: first, a previous reduction method for stability testing is recalled, and then, the proposed method is introduced and tested for several examples before concluding.

2 Phase stability testing using a reduction method

In conventional (in the compositional space) phase stability testing of an nc-component mixture of composition z_i, Michelsen’s modified TPD function [1]

is minimized with respect to formal mole numbers of the trial phase (Y_i), where φ_i is the fugacity coefficient.

A mixture is thermodynamically stable if D* is nonnegative for all feasible Y. The common procedure is to check the sign of D* at all calculated stationary points. In reduction methods, the nonideal part of the TPD function is expressed in terms of reduction parameters.

The M reduction parameters are defined as [10]

where q are the elements of the reduction matrix [10].

There are several procedures to achieve the reduction (to express the energy parameter A in the cubic EoS in terms of the reduction parameters of the form given by equation (2)) [9,10,11,12]. Here, the spectral decomposition of the matrix K with elements (1−kij);i,j=1,nc [11] is used (where k_ij is the BIP between components i and j) and A is reduced from a quadratic to a linear form:

and

where m is the rank of K (which is rank-deficient and usually m ≪ nc), λα;α=1,m are the eigenvalues of the matrix K and qαi;α=1,m;i=1,nc are the corresponding eigenvectors. It is worth noting that the eigenvalue-eigenvector problem has many applications in theoretical physics, engineering and chemistry (for instance, such eigenproblem has been translated for the alignment of molecules).

The last reduction parameter is the co-volume in the EoS (with M = m + 1 and qMi=Bi;i=1,nc).

Since the EoS parameters, A and B, as well as ∂A/∂xi are functions of reduction parameters, the fugacity coefficient is also a function of reduction parameters at P and T specifications, that is, φi=φi(P,T,Q), unlike in conventional methods where it is a function of phase composition, φi=φi(P,T,n).

The reduction parameters were first used as independent variables in phase stability testing by Nichita et al. (2006) [13] (this method is denoted here by NBH). The M error equations are

in M unknowns, Q=(Q1,Q2,...,QM)T, and the trial phase mole numbers, Y_i, are calculated from

where the index F stands for feed.

The Newton iteration equation is

where J(R) is the Jacobian matrix and F is the error vector given by equation (6),

From equation (7)

The solution of the linear system of equation (8) corresponds to a stationary point of D* [9]. The variables Q are bounded by mini=1,ncqαi≤Qα≤mini=1,ncqαi;α=1,M. Calculations start directly with Newton iterations; if a variable hits its bounds or if the TPD function increases between two subsequent Newton iterations, either a line search is preformed [13] or iterations are reverted to SSI [5]. SSI iterations are based on updating Y from equation (7) and Q from equations (2) and (5) using the normalised mole fractions.

Extensive testing showed [5,13] that the NBL reduction method is very fast in the two-phase region (convergence starting directly with Newton iterations is unproblematic even very close to critical points, where 10 iterations are rarely exceeded). In the single-phase region, the number of iterations exhibits a peak strictly near the singularity and an erratic behavior in a domain (that becomes wider with the distance to the critical point) outside the STLL, where convergence may be slow and even problematic.

3 Proposed method

Let us define the vector Q¯=(Q¯1,…,Q¯M,S)T of M + 1 modified reduction parameters, with elements

where Q¯α=SQα, as Yi=Sxi, and

with qM+1,i=1;i=1,nc.

It was shown that a solution of equation (11) corresponds to a stationary point of the modified TPD function D* and honors the constraints in its minimization with respect to a specific set of variables, i.e., Y and Q¯=(Q¯1,…,Q¯M,S)T [8].

The error equations are

The linear system of M + 1 equations is

The elements of the Jacobian matrix are

where

and

The partial derivatives of the fugacity coefficients with respect to the reduction parameters can be found in ref. [8] for a general form of two-parameter cubic EoS, containing Peng-Robinson [14] EoS and the Soave [15] EoS; these derivatives have simpler forms than those required in conventional methods.

The block structure of the Jacobian matrix (the prime on the symbol of an array means here that the array contains only its first M elements) is

Using simple linear algebra operations, it can be shown that an equivalent problem consists in first solving the linear system of M equations

for ΔQ¯α;α=1,M, where

and then updating Q¯M+1 (S) from

The Jacobian matrix J′ differs from its block counterpart in J by a rank one matrix (thus no additional summations are required), that is, the Schur complement of JM+1,M+1 is J; moreover, the condition number of J′ is not affected by the bordering vectors in equation (17). By this way, the dimensionality of the linear system resolution is not increased by adding an additional independent variable.

The algorithm for the proposed reduction method is as follows:

1. Provide initial estimates of the equilibrium constants (Wilson’s relation for ideal K-values is used [1])

2. Calculate initial estimates of formal mole numbers (Y_i)

3. Calculate initial values of the iteration variables from equations (11a) and (11b);

4. Calculate the fugacity coefficients as functions of the reduced variables

5. Update Y_i from equation (7)

6. Evaluate the error function (the Euclidean norm of the error vector in equations (12a) and (12b)

7. Convergence test; if the error function is less than the required tolerance, the iterative process is stopped.

8. Calculate the TPD function; if the TPD function is increasing between two subsequent iterations, an SSI iteration is performed (go to step iii)

9. Calculate the required partial derivatives and assemble the Jacobian matrix

10. Solve the linear system of equations (the Newton method)

11. Update the iteration variables and go to step iv

The iterative sequence above is performed twice, for vapor-like and liquid-like initialisations [1], as described in the next section.

The proposed method was coded in FORTRAN, according to the algorithm presented above, by modifications of the original code for the NBL reduction method for stability testing. The latter was previously implemented and thoroughly tested in several academic and industry research compositional simulators.

1. Ethical approval: The conducted research is not related to either human or animal use.

4 Results and discussion

Phase stability testing calculations are performed for two hydrocarbon mixtures. The PR EoS is used in all calculations. Calculations start directly with Newton iterations, with a possible switch-back to SSI when reduction parameters are going out of their bounds. The iterative process is stopped if the Euclidean norm

with ε=10−10 in all calculations.

Michelsen’s two-sided initialization is used, in which ideal equilibrium constants calculated with Wilson’s relation [1] are used to generate the initial guess. If the stability of a vapor-like mixture is investigated (the trial phase is liquid), initialization is denoted as type V and YiV(0)=zi/Ki(0); for a liquid-like mixture (the trial phase is vapor), the initialization is denoted as type L and YiL(0)=ziKi(0).

The first mixture is a six-component synthetic gas-condensate (Yarborough [16]), containing normal alkanes ranging from C₁ to nC₁₀, denoted as Y8. Composition is given in ref. [16], component properties are taken from Reid et al. [17], and all BIPs are set to zero. This test problem is a benchmark problem for phase equilibrium calculation, used in many studies [5,6,7,8,12,13,18], and covers the important particular case when there are only two iteration variables (the parameters of the EoS) for phase stability testing. There is only one nonzero eigenvalue of the matrix K, λ1=nc=6. The critical point calculated with the PR EoS using a reduction method [18] is T_c = 293.78 K and P_c = 210.67 bar. The second mixture (Metcalfe and Yarborough [19]), denoted as MY10, is a ten-component synthetic oil containing normal alkanes ranging from C₁ to nC₁₄; feed composition, component properties and nonzero BIPs are taken from ref. [5]. The matrix K is of rank 3; the three nonzero eigenvalues are λ1=9.957353, λ2=0.070650 and λ3=−0.028003. The calculated critical point using the PR EoS [18] is T_c = 572.23 K and P_c = 79.94 bar.

A large number of phase stability testing calculations are carried out with the proposed and NBH methods, by spanning the pressure–temperature plane using refined grids, with increments of ΔP = 1 bar and ΔT = 1 K. The temperature and pressure windows in the P−T plane are defined by [250 K, 600 K] and [1 bar, 300 bar], respectively, for the Y8 mixture (giving 105,350 test points) and by [250 K, 700 K] and [1 bar, 250 bar], respectively, for the MY10 mixture (giving 1,12,750 test points). The number of iterations for the NBH and proposed methods within the selected P−T windows are plotted in Figures 1 and 2 for the Y8 mixture with a type V initialization, Figures 3 and 4 for the MY10 mixture with a type L initialization and in Figures 5 and 6 for the MY10 mixture with a type V initialization. An increase in iteration numbers can be clearly identified in the vicinity of the STLL and in a certain domain in the single-phase region (in Figures 1–4) and in the vicinity of the spinodal curve in Figures 5 and 6. The explanation of this typical convergence behavior is given in ref. [6,7].

Figure 1

Number of iterations for the Y8 mixture, NBH method with type V initialization.

Figure 2

Number of iterations for the Y8 mixture, proposed method with type V initialization.

Figure 3

Number of iterations for the MY10 mixture, NBH method with type L initialization.

Figure 4

Number of iterations for the MY10 mixture, proposed method with type L initialization.

Figure 5

Number of iterations for the MY10 mixture, NBH method with type V initialization.

Figure 6

Number of iterations for the MY10 mixture, proposed method with type V initialization.

It appears from Figures 1–6 that the proposed method leads to a systematic decrease in the number of iterations in the single-phase region, mainly in the highly difficult P−T domains. In the two-phase region, the proposed method is slightly slower than NBH. Taking into account the above considerations, the recommended procedure is to start iterations with the proposed method and switch to NBH when the TPD function D* < 0 becomes negative (both the proposed and NBH [6,10,11,16] methods are highly robust in the two-phase region). We aimed in this paper to highlight the improvement of the original NBL method by using the proposed reduction method.

Finally, it is important to note that the results of the reduction method are identical to those obtained by any conventional methods. Differences with experimental values are those specific of the thermodynamic model (EoS) used. If the spectral expansion in equation (3) is truncated (by neglecting some small eigenvalues, as for instance in refs. [5,11]), noticeable differences may appear (see a detailed discussion in Michelsen et al. [20]).

As the proposed algorithm using the modified reduction parameters (Q¯) performs better than the original one using the reduction parameters (Q), the immediate practical consequence of this work is to extend the use of the new set of variables to more complex algorithms in a TPD function minimization framework.

5 Conclusions

A reduction method for mixture phase stability testing is proposed, using a new set of error functions and independent variables (the modified reduction parameters, including the sum of formal mole number) for Newton iterations. Numerical experiments carried out by spanning the pressure–temperature plane show the robustness of the proposed method, with no failure and a significantly improved convergence behavior, mainly in the single-phase region for pressures outside the STLL (where it is faster than the previous formulation), recommending the proposed method for compositional simulators, in which the problematic domain is crossed a very large number of times during simulations. A general form of the two-parameter cubic EoS is used, but the reduction method is applicable to any EoS in which all mixing rules are decomposable to linear forms.