Next Article in Journal
Entropy and Geometric Objects
Previous Article in Journal
Atom-Diffraction from Surfaces with Defects: A Fermatian, Newtonian and Bohmian Joint View
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Hybrid Newton–Successive Substitution Method for Multiphase Rachford-Rice Equations

1
School of Energy Resource, China University of Geosciences, Beijing 100083, China
2
Petroleum Engineering, Colorado School of Mines, Golden, CO 80401, USA
*
Authors to whom correspondence should be addressed.
Entropy 2018, 20(6), 452; https://doi.org/10.3390/e20060452
Submission received: 4 May 2018 / Revised: 2 June 2018 / Accepted: 6 June 2018 / Published: 9 June 2018
(This article belongs to the Section Thermodynamics)

Abstract

:
In multiphase (≥3) equilibrium calculations, when the Newton method is used to solve the material balance (Rachford-Rice) equations, poorly conditioned Jacobian can lead to false convergence. We present a robust successive substitution method that solves the multiphase Rachford-Rice equations sequentially using the method of bi-section while considering the monotonicity of the equations and the locations of singular hyperplanes. Although this method is slower than Newton solution, as it does not rely on Jacobians that can become poorly conditioned, it can be inserted into Newton iterations upon the detection of a poorly conditioned Jacobian. Testing shows that embedded successive substitution steps effectively improved the robustness. The benefit of the Newton method in the speed of convergence is maintained.

1. Introduction

Petroleum reservoir fluids are multi-component mixtures primarily made of hydrocarbons, held subsurface at high-temperature and high-pressure conditions [1,2]. Although under most situations the phase behavior of these fluids can be adequately described by vapor-liquid two-phase equilibria, multiphase phenomena (number of phases ≥3) have also been widely observed. Some of them, such as wax precipitation [3,4,5,6], hydrate formation [7,8,9,10], and asphaltene precipitation [11,12,13,14] are of significant importance to oil and gas production and transportation. In enhanced oil recovery, mixtures of reservoir oil and injected gas, such as CO2, can also exhibit complex phase behaviors at low pressures [15,16,17]. For these reasons, research on multiphase equilibrium calculations is still very active.
Multiphase (≥3) equilibrium calculation includes phase-stability [18] and phase-split [19] calculation steps. These two steps are usually carried out in series, except for the class of methods that was started by Gupta et al. [20], where they are carried out simultaneously. In the phase-split calculation step, mass balance equations and fugacity equations are solved. These two sets of equations can be solved sequentially using the method of successive substitution: the mass balance equations are first solved for a given set of equilibrium ratios; then, the equilibrium compositions of the phases that were obtained from the solution of the mass balance equations are used to update the equilibrium ratios through each phase’s equation of state or activity correlation [21,22,23]. These two sets of equations can also be solved simultaneously by the Newton method [24]. Although the Newton method converges more rapidly than the method of successive substitution, its accuracy is often influenced by the initial guess. For this reason, Michelson [19] recommended to begin phase-split calculations with successive substitutions and then switch to the Newton method. Successive substitution method and Newton method for phase-split calculations are abbreviated as SS-PS and NM-PS in this study.
This study focuses on the solution of the mass balance equations, the first step in SS-PS. Mass balance equations can be either solved by the Newton method [21,22,23,25,26] or by formulating the equations into a problem of minimization [27,28,29]. They are usually expressed in the Rachford-Rice form [30], which, for the case of two-phase equilibria, is easy to solve owing to its monotonicity and bounds of solutions explicitly indicated by the poles [31]. Following the abbreviations of SS-PS and NM-PS, we will abbreviate Newton method for the Rachford-Rice equations as NM-RR, and the minimization method for the Rachford-Rice equation as MM-RR. For multiphase (≥3) Rachford-Rice equations, the Newton method that was presented in Leibovici and Neoschil [26] is a widely-used NM-RR. Their method uses the locations of singular hyperplanes of the Rachford-Rice equations to confine the Newton iterations. The result is a robust convergence toward a unique solution in the multi-dimensional hypervolume, including the “negative flash” [32] region. It was found, however, by Okuno and coworkers [29], that Leibovici and Neoschil’s method would occasionally fail to converge to the correct answer, because the condition of the Jacobian becomes too large when the Newton steps are carried out very close to the singular hyperplanes. Okuno et al.’s method, which is based on the principle of minimization, uses the condition of non-negative phase compositions to further constrain the solutions and can overcome the problem of poorly conditioned Jacobians. Their method also converges with a fewer number of iterations as compared to the Leibovici and Neoschil’s method. The method of Okuno et al. is one example of the latest MM-RR.
In this study, we present a successive-substitution method to solve the multiphase Rachford-Rice equations that we will mention henceforth as SS-RR. Considerations of monotonicity and the locations of singular hyperplanes help this method to achieve robust convergence. The speed of convergence is slow when compared to NM-RR. However, as it does not use Jacobian or any derivatives, it can achieve convergence where NM-RR of Leibovici and Neoschil [26] cannot, due to poorly conditioned Jacobians. A hybrid code that combines the benefits of SS-RR and NM-RR is also presented. Section 2 is a presentation of the SS-RR method. In Section 3, we present the hybrid Newton-successive substitution method (NSS-RR) and two examples for which NM-RR could not converge to the correct answers, whereas both SS-RR and NSS-RR did. In Section 4, we present comparisons to show that the results from our calculations agree with those that were reported in the literature, including a complete phase-split calculation for an oil-gas-water system with equation of state and Henry’s law.

2. Successive Substitution Method for Multiphase Rachford-Rice Equations

In this study, we use n ˜ j to denote the mole fraction of phase j in a multiphase, multicomponent mixture, x i j to denote the mole fraction of component i in phase j, and z i to denote the mole fraction of component i in the mixture. N C is the total number of components and N P is the total number of phases. The equilibrium ratio of component i between phase α and phase β is defined as K i α β = x i α / x i β . The multiphase Rachford-Rice equations that are derived from mass balances are
i = 1 N C ( 1 K i j 1 ) z i 1 k = 2 N P n ˜ k ( 1 K i k 1 ) = 0 j = 2 ,   ,   N P
Here, selection of reference phase “1” is arbitrary, but it should contain every component of the system. While the selection of reference phase does not affect the results, it does affect convergence path and speed. We recommend to select reference phase to avoid extreme values of equilibrium ratios.
For a given set of equilibrium ratios K i j 1 ( i = 1 ,   2 ,   ,   N C , j = 2 ,   ,   N P ) and mole fractions z i , Equation (1) is solved to give n ˜ 2 through n ˜ N P . Then, the compositions of phases are obtained using
x i j = z i K i j 1 1 + k = 2 N P n ˜ k ( K i k 1 1 ) i = 1 ,   2 ,   ,   N C j = 1 ,   2 ,   ,   N P
In this study, we introduce
ξ i j = K i j 1 1 j = 2 ,   ,   N P
to simply Equation (1) into a set of fractional equations
F j ( n ˜ 2 , , n ˜ N P ) = i = 1 N C ξ i j z i 1 + k = 2 N P ξ i k n ˜ k = 0 j = 2 ,   ,   N P
The values of ξ i j vary between −1 and ∞.
Although Equation (4) appears simple, formulating a robust numerical solution is rather complex, because of the presence of N C hyperplanes
1 + k = 2 N P ξ i k n ˜ k = 0 i = 1 ,   2 ,   ,   N C
on which Equation (4) becomes singular. Precarious application of the Newton method, without proper consideration of the hyperplanes, can lead to non-convergence or convergence toward non-physical solutions. Leibovici and Neoschil [26] recognized that the physical constraint on the mole fractions of the phases, 0 n ˜ j 1 and j = 2 N n ˜ j < 1 , defines a hypertetrahedron in the [ n ˜ 2 ,   ,   n ˜ N P ] space, and this hypertetrahedron cannot be dissected by any of the hyperplanes. This property ensures that F j in Equation (4) are continuously differentiable in the hypervolume
1 + ξ i j n ˜ j > 0
that encloses the hypertetrahedron of the physically admissible solutions and the immediately adjacent negative flash region. By starting the initial guess within the hypervolume defined by Equation (6) and relaxing the Newton steps such that they do not cross the singular hyperplanes, Leibovici and Neoschil [26] were able to ensure that their Newton method does not generate a solution outside of Equation (6). However, the condition of F j being continuously differentiable does not guarantee that the Jacobian that is needed by the Newton method is always well conditioned. When Newton iteration approaches the boundaries of Equation (6), the Jacobian can become nearly singular, causing the Newton method to fail [29].
The abovementioned limitation of the Newton method motivated us to seek a method that does not require the computation of Jacobian or any derivatives of F j . Our method began with the recognition that since
F j n ˜ j = i = 1 N C ( ξ i j ) 2 z i [ 1 + k = 2 N P ξ i k n ˜ k ] 2 < 0
F j must be monotonically decreasing in the direction of increasing n ˜ j . Because of this property, when n ˜ j is increased to approach a singular hyperplane while all the other variables are kept constant, F j . When n ˜ j is decreased to approach a singular hyperplane while all the other variables are kept constant, F j + . Additionally, when n ˜ j + while all other variables are kept finite, F j decreases and approaches zero; when n ˜ j while all other variables are kept finite, F j increases and approaches zero. Figure 1 shows, qualitatively for a three-phase system, the monotonicity of F 2 and F 3 to n ˜ 2 and n ˜ 3 .
The monotonicity of F j to n ˜ j , but not the other variables motivated us to design a successive-substitution method to solve Equation (4), as follows. First, an initial guess is made within the hypertetrahedron of physically admissible solutions, n ˜ 0 2 = n ˜ 0 3 =     = n ˜ 0 N P = 1 / N P . Here, the subscripts to n ˜ j denote the number of iterative substitutions. We begin the first iteration by holding the values of n ˜ 0 3 through n ˜ 0 N P constant and seek n ˜ 1 2 to satisfy F 2 ( n ˜ 1 2 , n ˜ 0 3 ,   , n ˜ 0 N P ) = 0 . Once n ˜ 1 2 is identified, we hold n ˜ 1 2 and n ˜ 0 4 through n ˜ 0 N P constant and seek n ˜ 1 3 to satisfy F 3 ( n ˜ 1 2 , n ˜ 1 3 , n ˜ 0 4 ,   , n ˜ 0 N P ) = 0 . This process continues till F N P = 0 is solved and subscripts of all n ˜ j are updated to level “1”, which finishes the first iteration. The second iteration begins by holding the values of n ˜ 1 3 through n ˜ 1 N P constant and seek n ˜ 2 2 to satisfy F 2 ( n ˜ 2 2 , n ˜ 1 3 ,   , n ˜ 1 N P ) = 0 , and ends when all the subscripts of n ˜ j are updated to level “2”. Iteration would then begin at level “3”, and continue till all n ˜ j are converged.
With this successive-substitution strategy, solution of Equation (4) is reduced to successive solutions of single-variable equation F j ( n ˜ j ) = 0 . Recognizing that F j ( n ˜ j ) is monotonic between its poles, we used the following method to solve F j ( n ˜ j ) = 0 :
  • Evaluate F j ( n ˜ m 2 ,   ,   n ˜ m j 1 ,   [ n ˜ m 1 j * ] ,   n ˜ m 1 j + 1 ,   ,   n ˜ m 1 N P ) . Here, m refers to the level of iteration and n ˜ j * in the square bracket is the variable that needs to be updated to level “m” such that F j ( n ˜ m 2 ,   ,   n ˜ m j 1 ,   n ˜ m j ,   n ˜ m 1 j + 1 ,   ,   n ˜ m 1 N P ) = 0 .
  • Determine the direction along which to vary n ˜ j * . If F j ( n ˜ m 2 ,   ,   n ˜ m j 1 ,   [ n ˜ m 1 j * ] ,   n ˜ m 1 j + 1 ,   ,   n ˜ m 1 N P ) > 0 , n ˜ j * should be increased. If, however, F j ( n ˜ m 2 ,   ,   n ˜ m j 1 ,   [ n ˜ m 1 j * ] ,   n ˜ m 1 j + 1 ,   ,   n ˜ m 1 N P ) < 0 , n ˜ j * should be decreased.
  • Check whether a solution exists along the direction of increasing (or decreasing) n ˜ j * . This is achieved by performing a line search to see whether increasing (or decreasing) n ˜ j * would lead to an intersection with a singular hyperplane. If n ˜ j * needs to be increased and increasing n ˜ j * generates intersections with singular hyperplanes, one solution must exist between the current n ˜ j * and the nearest intersection, because, as n ˜ j * approaches the intersection, F j ( n ˜ j * ) approaches . If n ˜ j * needs to be increased, but there is no intersection with singular hyperplanes along the direction of increasing n ˜ j * , there is no solution to F j ( n ˜ m 2 ,   ,   n ˜ m j 1 ,   n ˜ m j ,   n ˜ m 1 j + 1 ,   ,   n ˜ m 1 N P ) = 0 and the calculation stops. Similar criteria apply when n ˜ j * needs to be decreased.
  • When solution exists, vary n ˜ j * to find n ˜ m j * . Because there is one and only one solution between n ˜ m 1 j * and the nearest intersection with singular hyperplanes, a bi-section method is used. Testing begins at the midpoint between n ˜ m 1 j * and the nearest intersection. The interval that contains the solution is continuously halved until F j evaluated at the midpoint of the interval becomes zero within a prescribed precision.
Successive substitution and bi-section methods that are presented above are very robust. As the initial guess is within the hypervolume defined by Equation (6) and solution search is bounded by the singular hyperplanes, this algorithm will never generate a solution outside the hypervolume. A flow chart for this algorithm is provided in Appendix B.
At the end of this section, we present an example that involves gas, oil, and water phases. The system that is considered here only consists of three components: methane, n-butane, and water. The pressure is 13.8 MPa and the temperature is 93.3 °C. Methane, n-butane, and moisture are present in the gas phase (g = 1); methane, n-butane, and dissolved water are present in the oil phase (o = 2); and, dissolved methane, dissolved n-butane, and water are present in the water phase (w = 3). The equilibrium ratios needed by the calculation were obtained from solubility data of pure hydrocarbons in pure water [33], published equilibrium ratios for hydrocarbons in a retrograde gas [34], moisture content in natural gas [35], and water solubility data in hydrocarbons [36]. At the pressure and temperature of interest, the mole fractions of methane and n-butane in the water phase are: x C 1 w = 1.6 × 10 3 , x n C 4 w = 1.5 × 10 4 . Assuming that the equilibrium ratios for methane and n-butane, between gas and oil phases, are not affected by whether these phases contain water or not, we obtained K C 1 g o = 2.181 , K n C 4 g o = 0.350 . Data on moisture content in natural gas suggests that at the pressure and temperature of interest, the mole fraction of water in the gas phase should be x H 2 O g = 2.894 × 10 2 . Correlation in Hibbard and Schalla [36] suggests that mole fraction of water in the oil phase can be set to x H 2 O o = 0.004 . The calculated compositions of the three phases and the equilibrium ratios are listed in Table 1, together with their corresponding ξ i j .
This example was solved using a mixture composition of z C 1 = 0.6 , z n C 4 = 0.35 and z H 2 O = 0.05 . The initial guess was ( n ˜ g ,   n ˜ o ,   n ˜ w ) = ( 1 / 3 ,   1 / 3 ,   1 / 3 ) . The case converged within five iterations with convergence criterion [ Δ n ˜ o , Δ n ˜ w ] 10 7 .
Figure 2 shows the path of convergence from the initial guess. Solution for n ˜ g , n ˜ o , and n ˜ w is (0.6725, 0.2981, 0.0294).
We note that this SS-RR method can be used to obtain, for three-phase equilibrium calculations, the lines on which F 2 = 0 and F 3 = 0 . These lines are of special interest, because their intersection is the needed solution. Such lines have been used by Haugen et al. [37] and Li and Firoozabadi [24] in order to explain their area-based bi-section method to solve F 2 and F 3 . The dotted lines in Figure 2 mark the locations where F 2 = 0 and F 3 = 0 , respectively. They were constructed from the paths of convergence of 1780 SS-RR calculations, each with a different initial guess within the area that is defined by Equation (6). These calculations all converged to the same solution, which indicates that initial guess is not important as long as it is bounded by Equation (6).

3. Hybrid Newton–Successive Substitution Method

The convergence speed of the above SS-RR method is slow when compared to the NM-RR method of Leibovici and Neochil [26]. However, as SS-RR is robust and it does not require the Jacobian or any derivative of F j , it can be integrated into NM-RR to handle situations with poorly conditioned Jacobians.
In hybrid Newton-Successive substitution (NSS-RR) method, the Newton step is identical to that in Leibovici and Neochil [26]. Prior to each Newton step, however, the condition number of the Jacobian matrix, κ , is evaluated. A Newton iteration is carried out if κ is less than a prescribed value and a successive-substitution iteration is carried out if otherwise. A flow chart for NSS-RR is provided in Appendix B.
When starting a success-substitution iteration, the result from the previous Newton step can be directly used. When starting a Newton step after a successive-substitution iteration, however, care should be taken because at the end of successive substitution F N P = 0 . Assume that [ n ˜ m 1 2 ,   ,   n ˜ m 1 N P ] is the result of the m 1 -th iteration and [ n ˜ m 2 ,   ,   n ˜ m N P ] is the result of the m -th iteration, which is a successive-substitution. If the m + 1 -th iteration is a Newton step, we recommend not to use [ n ˜ m 2 ,   ,   n ˜ m N P ] to start the Newton step, but a combination of [ n ˜ m 1 2 ,   ,   n ˜ m 1 N P ] and [ n ˜ m 2 ,   ,   n ˜ m N P ] . We specifically used
n ˜ i * = i 1 N P 1 n ˜ m 1 i + N P i N P 1 n ˜ m i
to start the Newton step after successive-substitution. Geometrically, n ˜ i * is the average of the points along the path of successive-substitution from [ n ˜ m 2 ,   n ˜ m 1 3 ,   ,   n ˜ m 1 N P ] to [ n ˜ m 2 , n ˜ m 3 , , n ˜ m N P ] . As each of these points makes one and only one F j zero, n ˜ j * calculated from Equation (8) is guaranteed to separate from surfaces F j = 0 , thus making it a good choice in our opinion to start the Newton step.
In what follows, we present two examples for which NM-RR of Leibovici and Neoschil [26] did not converge to the correct answers, whereas both SS-RR and NSS-RR were successful. The first example is a fifteen-component, three-phase mixture. The second example is a twenty-component, five-phase mixture. In NM-RR and NSS-RR, the relaxation parameter that regulates the Newton step when it intersects with the singular hyperplanes was set to 0.5. The maximum condition number in NSS-RR was set to 1010. In all of the methods, the convergence criterion set on the difference between vectors n ˜ m 1 i and n ˜ m i was 10−7. Note that if the convergence criterion is tightened, then NM-RR can find the correct solutions for these examples. We are only using these examples to compare NM-RR and NSS-RR under the same convergence criterion.
Table 2 gives the composition of the fifteen-component, three-phase mixture, the equilibrium ratios, and the equilibrium composition of the three phases. For this example, both NSS-RR and SS-RR converged to the correct root (−0.01686263294, −1.1254155641). NM-RR of Leibovici and Neoschil [26] however converged to a wrong root (−0.04078420653, −1.1004615900). As shown in Figure 3, many Newton steps were carried out very close to a singular line of the mixture. At the 21st step, a very large condition number of 5.0977 × 1011 was encountered. The hybrid method carried out a successive-substitution step at this location and the subsequent Newton steps converged to the correct solution. The NM-RR method, on the other hand, lost accuracy at this location and it converged to a wrong solution. Figure 4 shows the variations in the condition number of the Jacobian during iterations. NM-RR finished in 29 iteration steps and NSS-RR finished in 28 iteration steps. The 21st step in NSS-RR is the only successive-substitution step carried out.
Table 3 gives the composition of the twenty-component, five-phase mixture, the equilibrium ratios, and the equilibrium compositions of the five phases. For this example, both NSS-RR and SS-RR converged to the correct root (−0.00538660799, −0.00373696250, −0.00496311432, −0.00415370309). NM-RR of Leibovici and Neoschil [26], however, converged to a wrong root (−0.00287415017, −0.00392609623, −0.00798417906, −0.00350187286). Figure 5 shows that, at the 38th step, the condition number of the Jacobian reached 3.52 × 1011. Upon detecting this large condition number, NSS-RR performed a single successive-substitution step, which helped the following Newton steps to converge to the correct solution using a total of 54 iterations. NM-RR of Leibovici and Neoschil [26], however, lost accuracy after the 38th step and converged to a wrong solution.
As observed from the above examples, as SS-RR does not use any Jacobians nor derivatives, it is a useful method to switch to upon detection of a poorly conditioned Jacobian. In our hybrid NSS-RR, successive substitution is only activated on rare occasions, and hence it does not add significantly to the computational time. Checking the condition numbers of Jacobians, on the other hand, generated a non-negligible overhead to the algorithm. In average, the run time of NSS-RR is about 1.4 times of that of NM-RR.

4. Validations and A Three-Phase Equilibrium Calculation

After proving the robustness of the developed SS-RR and NSS-RR, in this section, we first compare our solutions of three-phase Rachford-Rice equations to those that are available in the literature. Then, we combine the developed solver with equation of state and Henry’s law to solve a complete oil-gas-water three-phase equilibrium.

4.1. Validations of Three-Phase Solutions

Three cases of three-phase equilibria from Nichita et al. [38] are used here to validate our solutions of three-phase Rachford-Rice equations. The first case is a ternary mixture containing CO2, CH4, and normal-hexadecane (nC16). The second case is a sour gas system with six components, the details of which can be found in Robinson et al. [39]. The third case is a quaternary mixture, the phase behavior of which was studied by Kohse and Heidemann [40]. Table 4, Table 5 and Table 6 list the compositions of the mixtures and the equilibrium compositions of the three phases.
Equilibrium ratios that were obtained from these tables were passed to our multiphase Rachford-Rice equation solvers. We note here that all solvers, NM-RR, SS-RR, and hybrid NSS-RR, converged on these cases and gave identical results. In Table 7, we compare our results to those that were reported in Nichita et al. [38]. This comparison shows that results from our algorithm are in very good agreement with previously reported solutions.
In a recent work by Okuno et al. [29], four three-phase examples were presented (Table 1 in reference [29]), including flash and negative flash near critical points. We note that our solvers can achieve convergence on all four examples. Although a quantitative comparison is not possible because the authors did not report the final mole fractions, our converged solutions do visually agree with that are those presented in the figures of reference [29].

4.2. Three-Phase Equilibrium of An Oil-Gas-Water System

In this section, we combine multiphase Rachford-Rice solvers with equation of state (EOS) and Henry’s law, and use the SS-PS approach to solve an oil-gas-water equilibrium. Water exists extensively in many hydrocarbon reservoirs. Water in reservoir may come from initial water traps, migration from nearby aquifers, or water injection for improved/enhanced oil recovery. The presence of water can affect the phase behavior of reservoir fluids, because some gases, in particular, CO2 and H2S, are soluble in the aqueous phase. Three-phase equilibrium calculations are needed to correctly describe the phase behaviors of such systems [41,42].
In this calculation, fugacities of components in the vapor and liquid phases were modeled using Peng-Robinson EOS following standard procedures. Fugacities of components in the water-rich phase were modeled using Henry’s law [43]. Henry’s law for a component sparingly soluble in the aqueous phase is
f i w = x i w H i , i w
where H i is Henry’s law constant of component i in the aqueous phase. Details of Henry’s law can be found in Appendix A. The fugacity of water in the aqueous phase can be obtained from the fugacity of the solutes using the Gibbs-Duhem equation, as in Prausnitz [44].
f H 2 O w = x H 2 O w Φ w s p v p w s ( p v p w s p v m R T d p )
Φ w s is the fugacity coefficient of pure water at the saturation vapor pressure p v p w s and v m is the molar volume of pure water. The fugacity coefficient of water Φ w s is calculated by the Chou equation, reported by Rowe and Chou [45], as follows:
Φ w s = { 0.9958 + 9.68330 × 10 5 T 6.175 × 10 7 T 2 3.08333 × 10 10 T 3 T > 90 F 1 T < 90 F
For a three-phase system with n c components, the following equations must be satisfied:
f i v = f i l i = 1 ,   n c
f i w = f i l i = 1 ,   n c
Li and Nghiem [46] defined the following sets of equilibrium ratios
K i v l = x i v x i l i = 1 ,   N c
K i v w = x i v x i w i = 1 ,   N c
that are consistent with our equilibrium ratios, as defined in Section 2. The three-phase Rachford-Rice equations for our vapor-liquid-water system, following Equation (1), are
i = 1 N C ( 1 K i v l ) z i 1 ( 1 K i v l ) n ˜ v ( 1 K i w l ) n ˜ w = 0
i = 1 N C ( 1 K i w l ) z i 1 ( 1 K i v l ) n ˜ v ( 1 K i w l ) n ˜ w = 0
The equilibrium calculation procedure began by passing initial guesses of equilibrium ratios to the Rachford-Rice equations. n ˜ l and n ˜ w solved from the Rachford-Rice equations were used to compute compositions of the phases: x i v , x i l , and x i w , following Equations (4) and (5). These compositions were then used to determine the fugacity coefficients Φ i v and Φ i l and the fugacity f i w , from which the equilibrium ratios were updated. The equilibrium ratios were again passed to the Rachford-Rice equations, until the results converge.
We applied the method that is presented above to a six-component mixture, the composition of which is shown in Table 8. Detailed description of this case can be referred to Li and Nghiem [46]. We conducted equilibrium calculation at 10 MPa and 100 °C. Our calculated results are in good agreement with results in the literature, as shown in Table 9.
To illustrate the difference between a three-phase equilibrium calculation and two-phase equilibrium calculations that discount the influence of the aqueous phase, we conducted a two-phase equilibrium calculation where water is removed from the mixture. The mole fractions of the non-aqueous components are reconstructed based on the relative fractions of them in the original mixture. The results of this two-phase calculation are also shown in Table 9. Comparison shows that although the two-phase approach neglecting the water component generated phase compositions that are similar to those from the three-phase calculation, the ratio between liquid and gas mole fractions, however, is noticeably higher. At equilibrium, the aqueous phase contains 14% of H2S and 2% of CO2 by mole, and is therefore highly corrosive. This corrosive nature can only be captured by a three-phase equilibrium calculation.

5. Conclusions

This paper presents a successive-substitution method to solve multiphase Rachford-Rice equations. In this algorithm, Rachford-Rice equations are solved sequentially, and the solution of each individual Rachford-Rice equation is achieved using the method of bi-section. The direction of bi-section is determined by the monotonicity of the equation, whereas the limit of bi-section is determined by the locations of the poles. These considerations ensure that the algorithm always reliably converges to either a physically admissible solution or a solution in the adjacent negative flash region.
In compositional reservoir simulations, the number of phase equilibrium calculations performed is very large. Accuracy and robustness without a significant penalty on efficiency are therefore very important. A hybrid method was developed in order to combine the advantages of this successive-substitution method and the Newton method. In this hybrid method, a successive substitution step replaces a Newton step when the local Jacobian becomes poorly conditioned. This hybrid method can effectively suppress the errors in the Newton method of Leibovici and Neochil [26] owing to poorly conditioned Jacobian. In terms of computational time, successive substitution steps add very little to the overall computational cost. The computational time of this hybrid method is 1.4 times of that of the Newton method of Leibovici and Neochil [26], and this difference is mainly due to the overhead needed to check the condition numbers.
We presented seven examples to show the characteristics and the accuracy of the algorithm. The first example shows the convergence behavior of the successive-substitution method. The second and the third examples present the advantages of the hybrid method over the Newton method of Leibovici and Neochil [26]. In the rest of the examples, we verified our algorithm using results from the literature. Lastly, we combined the developed multiphase Rachford-Rice solvers with equations of states for the oil and gas phases and Henry’s law for the water phase and completed a three-phase flash calculation for an oil-gas-water system with sour gases. Our result again is in very good agreement with that previously reported. It shows specifically that a two-phase equilibrium calculation neglecting the dissolution of sour gases in the water phase is not a good approximation of a true three-phase flash calculation.

Author Contributions

R.G. performed the investigation, wrote the original manuscript and made the revisions. X.Y. and Z.L. are supervisors of this study.

Funding

R.G. and Z.L. thank the financial support from the National Science and Technology Major Project under Grant NO. 2017ZX05009005. R.G. specifically thanks the China Scholarship Council for sponsoring R.G.’s visit to the Colorado School of Mines.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

A ( T ) H ( T ) functions used to compute the specific volume v
C cohesive energy density of water
F Rachford-Rice equations defined in Equation (4)
f i j fugacity of component i in phase j
h enthalpy departure
H i Henry’s law constant of component i
K i α β equilibrium ratio of component i between phase α and phase β
n ˜ j mole fraction of phase j
N c number of components
N P total number of phases
p pressure
R universal gas constant
T temperature
v specific volume
v m molar volume of pure water
w N a C l weight fraction of NaCl
x i j mole fraction of component i in phase j
z i mole fraction of component i in the multiphase mixture
κ condition number of Jacobian
ξ i j defined in Equation (3) as K i j 1 1
Φ i fugacity coefficient of component i
Superscripts
ggas
jphase index
lliquid
ooil
ssaturated
wwater
wswater at saturation
α , β phase index
Subscripts
ccritical point
i , j , k component index
vpvapor pressure
mlevel of iteration

Appendix A. Henry’s Law

Henry’s law constant can be solved from a differential equation proposed by Smith and Van Ness [47]
d ( ln H i ) = v m i R T d p + h i v h i R T 2 d T
In this equation, v m i is the partial molar volume of component i in the aqueous phase at infinite dilution, h i v is the enthalpy of component i in the gas phase, and h i is the enthalpy of component i in the aqueous phase at infinite dilution. For a given temperature T, integration of Equation (A1) from p 0 to p gives
ln H i = ln H i 0 + v m i ( p p i 0 ) R T
where H i 0 is Henry’s law constant at the reference pressure p i 0 . Equation (A2) can also be written as
ln H i = ln H i * + v m i p / ( R T )
where
ln H i * = ln H i 0 v m i p i 0 / ( R T )
In this study, Equation (A3) was used to correlate Henry’s law constant from solubility data. H i * is referred to as the reference Henry’s law constant. The molar volume at infinite dilution was calculated using the following correlation proposed by Lyckman et al. [48]
p c i v m i R T c i = 0.095 + 2.35 ( T p c i C T c i )
where C is the cohesive energy density of water given by:
C = ( h w 0 h w s + p v p w s v m w s R T ) v m w s
p v p w s is the saturated vapor pressure of water at temperature T, v m w s is the molar volume of water at p v p w s and T, and h w o h w s is the enthalpy departure of liquid water at p v p w s and T. The saturated vapor pressure of water p v p w s was calculated using the Buck equation [49]
p v p w s = 0.61121 exp ( ( 18.678 T 234.5 ) ( T 257.14 + T ) )
T is in C , and p v p w s is in p s i . Equation (A7) is applicable to liquid water: T > 0   C . When T < 0   C , for ice we used the following equation
p v p w s = 0.61115 exp ( ( 23.036 T 333.7 ) ( T 279.82 + T ) )
Enthalpy departure of liquid water ( h w 0 h w s ) were calculated from Yen-Alexander, reported in Poling et al. [50] in the form of:
h w 0 h w s = T c 7 4.5688 ( ln ( p v p w s 217.6 ) ) 0.333 1 + 0.004 ln ( p v p w s 217.6 )
The molar volume v m was estimated from the specific volume v , which, in turn, was from a correlation by Chou reported in Rowe and Chou [45].
v = A ( T ) p B ( T ) p 2 C ( T ) + w N a C l D ( T ) + x 2 E ( T ) x p F ( T ) x 2 p G ( T ) 1 2 x p 2 H ( T )
v is the specific volume in cm 3 / g , T is the temperature in K , p is the absolute pressure in kgf / cm 2 , and w N a C l is the weight fraction of NaCl in solution. The temperature functions are,
A ( T ) = 5.916365 0.01035794 T + 0.9270048 × 10 5 T 2 1127.522 / T + 100674.1 / T 2
B ( T ) = 0.5204914 × 10 2 0.10482101 × 10 4 T + 0.8328532 × 10 8 T 2 1.1702939 / T + 102.2783 / T 2
C ( T ) = 0.118547 × 10 7 0.6599143 × 10 10 T
D ( T ) = 2.5166 + 0.0111766 T 0.170552 × 10 4 T 2
E ( T ) = 2.84851 0.0154395 T + 0.223982 × 10 4 T 2
F ( T ) = 0.0014814 + 0.82969 × 10 5 T 0.12469 × 10 7 T 2
G ( T ) = 0.0027141 0.15391 × 10 4 T + 0.22655 × 10 7 T 2
H ( T ) = 0.62158 × 10 5 0.40075 × 10 8 T + 0.65972 × 10 11 T 2

Appendix B. Flow Charts

Figure A1. Flow chart of SS-RR.
Figure A1. Flow chart of SS-RR.
Entropy 20 00452 g0a1
Figure A2. Flow chart of NSS-RR.
Figure A2. Flow chart of NSS-RR.
Entropy 20 00452 g0a2

References

  1. Whitson, C.H.; Brulé, M.R. Phase Behavior; Society of Petroleum Engineers: Houston, TX, USA, 2000. [Google Scholar]
  2. Pedersen, K.S.; Christensen, P.L.; Shaikh, J.A. Phase Behavior of Petroleum Reservoir Fluids, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2014. [Google Scholar]
  3. Won, K.W. Thermodynamics for solid solution-liquid-vapor equilibria: Wax phase formation from heavy hydrocarbon mixtures. Fluid Phase Equilib. 1986, 30, 265–279. [Google Scholar] [CrossRef]
  4. Pedersen, K.S. Prediction of cloud point temperatures and amount of wax precipitation. SPE Prod. Facil. 1995, 10, 46–49. [Google Scholar] [CrossRef]
  5. Lira-Galeana, C.; Firoozabadi, A.; Prausnitz, J.M. Thermodynamics of wax precipitation in petroleum mixtures. AIChE J. 1996, 42, 239–248. [Google Scholar] [CrossRef] [Green Version]
  6. Coutinho, J.A. Predictive UNIQUAC: A new model for the description of multiphase solid-liquid equilibria in complex hydrocarbon mixtures. Ind. Eng. Chem. Res. 1998, 37, 4870–4875. [Google Scholar] [CrossRef]
  7. Bishnoi, P.R.; Gupta, A.K.; Englezos, P.; Kalogerakis, N. Multiphase equilibrium flash calculations for systems containing gas hydrates. Fluid Phase Equilib. 1989, 53, 97–104. [Google Scholar] [CrossRef]
  8. Ballard, A.L.; Sloan, E.D., Jr. The next generation of hydrate prediction: Part III. Gibbs energy minimization formalism. Fluid Phase Equilib. 2004, 218, 15–31. [Google Scholar] [CrossRef]
  9. Segtovich, I.S.V.; Barreto, A.G., Jr.; Tavares, F.W. Phase diagrams for hydrates beyond incipient condition-Complex behavior in methane/propane and carbon dioxide/iso-butane hydrates. Fluid Phase Equilib. 2016, 426, 75–82. [Google Scholar] [CrossRef]
  10. 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 Phsase Equilib. 2016, 414, 117–132. [Google Scholar] [CrossRef]
  11. Thomas, F.B.; Bennion, D.B.; Bennion, D.W.; Hunter, B.E. Experimental and theoretical studies of solids precipitation from reservoir fluid. J. Can. Pet. Technol. 1992, 31. [Google Scholar] [CrossRef]
  12. MacMillan, D.J.; Tackett, J.E., Jr.; Jessee, M.A.; Monger-McClure, T.G. A unified approach to asphaltene precipitation: Laboratory measurement and modeling. J. Pet. Technol. 1995, 47, 788–793. [Google Scholar] [CrossRef]
  13. Nghiem, L.X.; Coombe, D.A. Modelling asphaltene precipitation during primary depletion. SPE J. 1997, 2, 170–176. [Google Scholar] [CrossRef]
  14. Almehaideb, R.A. Asphaltene precipitation and deposition in the near wellbore region: A modeling approach. J. Pet. Sci. Eng. 2004, 42, 157–170. [Google Scholar] [CrossRef]
  15. Orr, F.M., Jr.; Yu, A.D.; Lien, C.L. Phase behavior of CO2 and crude oil in low-temperature reservoirs. Soc. Pet. Eng. J. 1981, 21, 480–492. [Google Scholar] [CrossRef]
  16. Henry, R.L.; Metcalfe, R.S. Multiple-phase generation during carbon dioxide flooding. Soc. Pet. Eng. J. 1983, 23, 595–601. [Google Scholar] [CrossRef]
  17. Turek, E.A.; Metcalfe, R.S.; Fishback, R.E. Phase behavior of several CO2/west Texas-Reservoir-Oil Systems. SPE Reserv. Eng. 1988, 3, 505–516. [Google Scholar] [CrossRef]
  18. Michelsen, M.L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1–19. [Google Scholar] [CrossRef]
  19. Michelsen, M.L. The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilib. 1982, 9, 21–40. [Google Scholar] [CrossRef]
  20. Gupta, A.K.; Bishnoi, P.R.; Kalogerakis, N. A method for the simultaneous phase equilibria and stability calculations for multiphase reacting and non-reacting systems. Fluid Phase Equilib. 1991, 63, 65–89. [Google Scholar] [CrossRef]
  21. Mehra, R.K.; Heidemann, R.A.; Aziz, K. Computation of multiphase equilibrium for compositional simulation. Soc. Pet. Eng. J. 1982, 22, 61–68. [Google Scholar] [CrossRef]
  22. Nghiem, L.X.; Li, Y.K. Computation of multiphase equilibrium phenomena with an equation of state. Fluid Phase Equilib. 1984, 17, 77–95. [Google Scholar] [CrossRef]
  23. Risnes, R. Equilibrium calculations for coexisting liquid phases. Soc. Pet. Eng. J. 1984, 24, 87–96. [Google Scholar] [CrossRef]
  24. Li, Z.; Firoozabadi, A. Initialization of phase fractions in Rachford–Rice equations for robust and efficient three-phase split calculation. Fluid Phase Equilib. 2012, 332, 21–27. [Google Scholar] [CrossRef]
  25. Nelson, P.A. Rapid phase determination in multiple-phase flash calculations. Comput. Chem. Eng. 1987, 11, 581–591. [Google Scholar] [CrossRef]
  26. Leibovici, C.F.; Neoschil, J. A solution of Rachford-Rice equations for multiphase systems. Fluid Phase Equilib. 1995, 112, 217–221. [Google Scholar] [CrossRef]
  27. Michelsen, M.L. Calculation of multiphase equilibrium. Comput. Chem. Eng. 1994, 18, 545–550. [Google Scholar] [CrossRef]
  28. Leibovici, C.F.; Nichita, D.V. A new look at multiphase Rachford–Rice equations for negative flashes. Fluid Phase Equilib. 2008, 267, 127–132. [Google Scholar] [CrossRef]
  29. Okuno, R.; Johns, R.; Sepehrnoori, K. A new algorithm for rachford-rice for multiphase compositional simulation. SPE J. 2010, 15, 313–325. [Google Scholar] [CrossRef]
  30. Rachford, H.H., Jr.; Rice, J.D. Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium. J. Pet. Technol. 1952, 4. [Google Scholar] [CrossRef]
  31. Leibovici, C.; Neoschil, J. A new look at the Rachford-Rice equation. Fluid Phase Equilib. 1992, 74, 303–308. [Google Scholar] [CrossRef]
  32. Whitson, C.H.; Michelsen, M.L. The negative flash. Fluid Phase Equilib. 1989, 53, 51–71. [Google Scholar] [CrossRef]
  33. Brooks, W.B.; Gibbs, G.B.; McKetta, J.J., Jr. Mutual solubilities of light hydrocarbonwater systems. Pet. Refin. 1951, 30, 118–120. [Google Scholar]
  34. Roland, C.H.; Smith, D.E.; Kaveler, H.H. Equilibrium constants for a gas-distillate System. Oil Gas J. 1941, 39, 128. [Google Scholar]
  35. Bukacek, R.F. Equilibrium Moisture Content of Natural Gases; Institute of Gas Technology: Des Plaines, IL, USA, 1955. [Google Scholar]
  36. Hibbard, R.R.; Schalla, R.L. Solubility of Water in Hydrocarbons; NACA-RM-E52D24; National Advisory Committee for Aeronautics: Washington, WA, USA, 10 July 1952.
  37. Haugen, K.B.; Firoozabadi, A.; Sun, L. Efficient and robust three-phase split computations. AIChE J. 2011, 57, 2555–2565. [Google Scholar] [CrossRef]
  38. Nichita, D.V.; Gomez, S.; Luna, E. Multiphase equilibria calculation by direct minimization of Gibbs free energy with a global optimization method. Comput. Chem. Eng. 2002, 26, 1703–1724. [Google Scholar] [CrossRef]
  39. Robinson, D.B.; Kalra, H.; Rempis, H. The Equilibrium Phase Properties of A Synthetic Sour Gas Mixture and a Simulated Natural Gas Mixture; Gas Processors Association: Tulsa, OK, USA, 1978. [Google Scholar]
  40. Kohse, B.F.; Heidemann, R.A. Tricritical lines and multiphase equilibria in quaternary mixtures. Fluid Phase Equilib. 1992, 75, 11–22. [Google Scholar] [CrossRef]
  41. Erbar, J.H.; Jagota, A.K.; Muthswamy, S.; Moshfeghian, M. Predicting Synthetic Gas and Natural Gas Thermodynamic Properties Using a Modified Soave-Redlich-Kwong Equation of State; Gas Processor Research Report, GPA RR-42; Gas Processors Association: Tulsa, OK, USA, 1980. [Google Scholar]
  42. Reshadi, P.; Nasrifar, K.; Moshfeghian, M. Evaluating the phase equilibria of liquid water+ natural gas mixtures using cubic equations of state with asymmetric mixing rules. Fluid Phase Equilib. 2011, 302, 179–189. [Google Scholar] [CrossRef]
  43. 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. [Google Scholar] [CrossRef]
  44. Prausnitz, J.M.; Lichtenthaler, R.N.; de Azevedo, E.G. Molecular Thermodynamics of Fluid-Phase Equilibria; Pearson Education: London, UK, 1998. [Google Scholar]
  45. Rowe, A.M., Jr.; Chou, J.C. Pressure-volume-temperature-concentration relation of aqueous sodium chloride solutions. J. Chem. Eng. Data 1970, 15, 61–66. [Google Scholar] [CrossRef]
  46. Li, Y.K.; Nghiem, L.X. Phase equilibria of oil, gas and water/brine mixtures from a cubic equation of state and Henry’s law. Can. J. Chem. Eng. 1986, 64, 486–496. [Google Scholar] [CrossRef]
  47. Smith, J.M.; Van Ness, H.C. Introduction to Chemical Engineering Thermodynamics, 3rd ed.; McGraw-Hill: Blacklick, OH, USA, 1975; pp. 241–244. [Google Scholar]
  48. Lyckman, E.W.; Eckert, C.A.; Prausnitz, J.M. Generalized reference fugacities for phase equilibrium thermodynamics. Chem. Eng. Sci. 1965, 20, 685–691. [Google Scholar] [CrossRef]
  49. Buck, A.L. New equations for computing vapor pressure and enhancement factor. J. Appl. Meteorol. 1981, 20, 1527–1532. [Google Scholar]
  50. Poling, B.E.; Prausnitz, J.M.; O’Connell, J.P. The Properties of Gases and Liquids, 5th ed.; McGraw-Hill: Blacklick, OH, USA, 2000. [Google Scholar]
Figure 1. (a) The blue triangle is the domain of physically admissible solutions, surrounded by four singular lines. The variation in F 2 along the horizontal dashed line A-B-C and that in F 3 along the vertical dashed line D-E-F are shown in (b,c), respectively. Both F 2 and F 3 decrease monotonically with increasing n ˜ 2 and n ˜ 3 . The monotonic trends are separated by poles at A, B, C and D, E, F.
Figure 1. (a) The blue triangle is the domain of physically admissible solutions, surrounded by four singular lines. The variation in F 2 along the horizontal dashed line A-B-C and that in F 3 along the vertical dashed line D-E-F are shown in (b,c), respectively. Both F 2 and F 3 decrease monotonically with increasing n ˜ 2 and n ˜ 3 . The monotonic trends are separated by poles at A, B, C and D, E, F.
Entropy 20 00452 g001
Figure 2. Red line ended with circle and cross marks the path of convergence for the three-component, three-phase example presented in Table 1. Circle is the initial guess and cross is the converged solution. The smaller shaded triangle is the domain of the physically permissible solutions. The larger shaded triangle is the domain defined by Equation (6) including the negative flash region. The dashed lines are the singular lines. The dotted lines inside the larger shaded triangle correspond to F 2 = 0 and F 3 = 0 , the intersection of which is the converged solution.
Figure 2. Red line ended with circle and cross marks the path of convergence for the three-component, three-phase example presented in Table 1. Circle is the initial guess and cross is the converged solution. The smaller shaded triangle is the domain of the physically permissible solutions. The larger shaded triangle is the domain defined by Equation (6) including the negative flash region. The dashed lines are the singular lines. The dotted lines inside the larger shaded triangle correspond to F 2 = 0 and F 3 = 0 , the intersection of which is the converged solution.
Entropy 20 00452 g002
Figure 3. Paths of convergence for the case presented in Table 2. Red line with dots is the path of SS-RR, blue line with dots is the path of NM-RR, and blue line with circles is the path of the hybrid method. The graph on the left shows the entirety of the paths and the graph on the right focuses on the final convergence.
Figure 3. Paths of convergence for the case presented in Table 2. Red line with dots is the path of SS-RR, blue line with dots is the path of NM-RR, and blue line with circles is the path of the hybrid method. The graph on the left shows the entirety of the paths and the graph on the right focuses on the final convergence.
Entropy 20 00452 g003
Figure 4. Variations in the condition number of Jacobian during NM-RR and NSS-RR iterations for the fifteen-component, three-phase example.
Figure 4. Variations in the condition number of Jacobian during NM-RR and NSS-RR iterations for the fifteen-component, three-phase example.
Entropy 20 00452 g004
Figure 5. Variations in the condition number of Jacobian during NM-RR and NSS-RR iterations for the twenty-component, five-phase example.
Figure 5. Variations in the condition number of Jacobian during NM-RR and NSS-RR iterations for the twenty-component, five-phase example.
Entropy 20 00452 g005
Table 1. Compositions of phases and equilibrium ratios.
Table 1. Compositions of phases and equilibrium ratios.
x i g x i o x i w K i g o K i g w ξ i o ξ i w
H2O0.028940.0040.998257.2350.029−0.86233.494
CH40.741430.339950.00162.181463.394−0.541−0.998
n-C4H100.229630.656050.000150.3501530.8671.857−0.999
Table 2. Composition and equilibrium ratios of the 15-component, 3-phase example.
Table 2. Composition and equilibrium ratios of the 15-component, 3-phase example.
Component z i K i 1 K i 2 x i 1 x i 2 x i 3
10.03155833298031.85287414286631.81156597622430.43667679408100.80910714054240.7910688227638
20.40712700766230.23140553083860.69549098601570.30031652088730.06949490393550.2088674332287
30.47519416717260.50417094443350.00010845017670.22271373742790.11228579533730.0000241533442
40.05458117115660.06354820838970.00126035117200.02550774486000.00162097148590.0000321487161
50.01157004468950.40789085062720.00134746944810.00542205984420.00221160860200.0000073060600
60.01891139557040.50662314810750.00000389293190.00886306530860.00449023404850.0000000345033
70.000045548480527.19016896435800.00352191331660.00002711514860.00073726547060.0000000954972
80.00064040143740.07650956934230.00001719238360.00029911760910.00002288535950.0000000051425
90.00036759736870.12849928328370.00000219653000.00017176573150.00002207177340.0000000003773
100.00000375048951.47955578722480.00016338404360.00000177148440.00000262101000.0000000002894
110.000000242884612.77698842934170.00160902285360.00000012617290.00000161211000.0000000002030
120.000000159440813.76668446722690.00075230461700.00000008350800.00000114962770.0000000000628
130.000000022858952.49955619490130.00007986824010.00000001818660.00000095478790.0000000000015
140.000000020254333.95392406721090.00000235160370.00000001290310.00000043811000.0000000000000
150.00000013755375.19791943338210.00001275744890.00000006694870.00000034799390.0000000000009
Table 3. Composition and equilibrium ratios of the 20-component, 5-phase example.
Table 3. Composition and equilibrium ratios of the 20-component, 5-phase example.
Component z i K i 1 K i 2 K i 3 K i 4 x i 1 x i 2 x i 3 x i 4 x i 5
10.38173995091402.37889143187140.18263462522182.13411484333780.71010732361420.38512819219760.91617815659100.07033774304440.82190779155680.2734823498098
20.07643364337310.83545374044022.06842866859201.90430183929436.04408598953890.07867964192240.06573320114060.16274322698690.14982978682790.4755465214053
30.13914877375700.11559384612542.84731834761620.01449452097990.43690411602930.13844399699430.01600327408560.39419413275930.00200667941900.0604867521264
40.06439922189520.00622628306252.13838603819280.04421689367810.99184888669950.06402299195860.00039862527040.13690587212760.00283089782840.0635011332973
50.14860260049510.00221565842480.79464161113260.07873371700420.77688845551860.14689259751700.00032546382120.11672697035430.01156540020290.1141191632121
60.04172124866530.01159514447652.16034343679410.05604949509960.21346117955370.04135156679220.00047947739130.08933358592060.00231773444030.0088269542238
70.12276935007670.00641674722550.15937920345960.07700424127530.02399489656880.12070192160390.00077451372060.01923737612130.00929455989360.0028962301245
80.02130878702390.00389463210180.03359176241380.00250502311280.02180594214170.02093219856550.00008152321250.00070314944100.00005243564120.0004564463108
90.00162703503090.01343664967200.72232584159190.10317431670400.17080861193880.00160418003540.00002155480510.00115874069410.00016551017900.0002740077651
100.00213074323060.00087340249972.61327064802390.00221309570420.09327274959550.00211368247830.00000184609560.00552362437980.00000467778160.0001971489765
110.00009178103050.010884487033324.40650053095080.19286907291871.00144148816360.00009956084880.00000108366880.00242993190850.00001920220860.0000997043646
120.00002298319300.030528838588125.84948987909190.05883930756724.09965906708580.00002541948340.00000077602730.00065708067890.00000149566480.0001042112156
130.00000347825510.018420675849210.47488595518600.35568523361810.10453828191990.00000356087680.00000006559380.00003729977800.00000126655130.0000003722479
140.00000011263671.955694412375657.64251280904231.748677793271829.05784702003480.00000016990950.00000033229110.00000979401170.00000029711700.0000049372049
150.00000023446340.28744670367821.04191876604361.888571945937313.70026993111250.00000024771140.00000007120380.00000025809510.00000046782070.0000033937126
160.00000000380641.535677537300653.551391118356597.73610710360556.64835339429090.00000001280930.00000001967100.00000068595660.00000125193260.0000000851609
170.00000001731260.75742722307867.69104012879616.007223802222918.77420855741800.00000001972670.00000001494150.00000015171880.00000011850270.0000003703530
180.00000002813660.00743777137576.76817274780284.05747619827245.27792810967420.00000002959110.00000000022010.00000020027750.00000012006510.0000001561795
190.00000000425890.000457402402928.139411565950935.15531735217789.05400327597300.00000000607550.00000000000280.00000017096060.00000021358560.0000000550075
200.00000000244530.08475613306131.648649403362531.96763170624802.51584408110750.00000000290230.00000000024600.00000000478490.00000009278110.0000000073018
Table 4. Composition of CO2–CH4–nC16 ternary mixture [38] and equilibrium compositions of the three phases (g: gas; l1 and l2: liquids).
Table 4. Composition of CO2–CH4–nC16 ternary mixture [38] and equilibrium compositions of the three phases (g: gas; l1 and l2: liquids).
Component z i x i g x i l 1 x i l 2
C10.050.0781120.0361810.038707
nC160.050.0000690.3402240.004609
CO20.900.9218190.6235950.956683
Table 5. Composition of the sour gas system mixture [38,39], and equilibrium compositions of the three phases (g: gas; l1 and l2: liquids).
Table 5. Composition of the sour gas system mixture [38,39], and equilibrium compositions of the three phases (g: gas; l1 and l2: liquids).
Component z i x i g x i l 1 x i l 2
C10.705920.7261290.7149730.065262
C20.068600.0063520.0721540.012393
C30.029670.0004660.0313330.003738
H2S0.105590.0085960.0974940.897742
CO20.019960.0049330.0206210.019166
N20.070260.2535240.0634250.001699
Table 6. Composition of the quaternary mixture [38,40] and equilibrium compositions of the three phases (g: gas; l1 and l2: liquids).
Table 6. Composition of the quaternary mixture [38,40] and equilibrium compositions of the three phases (g: gas; l1 and l2: liquids).
Component z i x i g x i l 1 x i l 2
C10.500000.9775990.3086460.068425
nC60.026270.0000010.2272520.007764
H2S0.416330.0079000.3634830.833419
CO20.057400.0144990.1006180.090391
Table 7. Mole fractions of the three phases (g: gas; l1 and l2: liquids) from this study and from Nichita et al. [38].
Table 7. Mole fractions of the three phases (g: gas; l1 and l2: liquids) from this study and from Nichita et al. [38].
The Ternary MixtureThe Sour Gas SystemThe Quaternary Mixture
This Study[38]DifferencesThis Study[38]DifferencesThis Study[38]Differences
n ˜ g 0.29080.29621.8%0.04170.04072.5%0.44840.44820.0%
n ˜ l 1 0.57000.56451.0%0.94340.94470.1%0.10090.10020.7%
n ˜ l 2 0.13920.13930.1%0.01490.01462.1%0.45070.45160.2%
Table 8. Composition of the six-component mixture.
Table 8. Composition of the six-component mixture.
ComponentComposition (%)
C130
nC515
nC1025
CO210
H2S
H2O
10
10
Table 9. Results for the case with three phases and six components in presence of water.
Table 9. Results for the case with three phases and six components in presence of water.
This StudyLi and Nghiem (1986)Two-Phase
x l % x g % x w % x l % x g % x w % x l % x g %
C122.86865.1740.00122.88465.2200.00123.28666.365
nC521.6764.936020.2164.603020.5854.687
nC1035.9190.663035.8890.662035.9590.664
CO29.14516.8730.0209.14816.8790.0219.23917.047
H2S10.90611.2200.14110.90611.2190.14110.92911.243
H2O0.9581.41699.7730.9581.41799.839//
n ˜ l n ˜ g n ˜ w n ˜ l n ˜ g n ˜ w n ˜ l n ˜ g
69.2021.759.0569.2621.709.0477.4522.55

Share and Cite

MDPI and ACS Style

Gao, R.; Yin, X.; Li, Z. Hybrid Newton–Successive Substitution Method for Multiphase Rachford-Rice Equations. Entropy 2018, 20, 452. https://doi.org/10.3390/e20060452

AMA Style

Gao R, Yin X, Li Z. Hybrid Newton–Successive Substitution Method for Multiphase Rachford-Rice Equations. Entropy. 2018; 20(6):452. https://doi.org/10.3390/e20060452

Chicago/Turabian Style

Gao, Ran, Xiaolong Yin, and Zhiping Li. 2018. "Hybrid Newton–Successive Substitution Method for Multiphase Rachford-Rice Equations" Entropy 20, no. 6: 452. https://doi.org/10.3390/e20060452

APA Style

Gao, R., Yin, X., & Li, Z. (2018). Hybrid Newton–Successive Substitution Method for Multiphase Rachford-Rice Equations. Entropy, 20(6), 452. https://doi.org/10.3390/e20060452

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop