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 CO
2, 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
to denote the mole fraction of phase
j in a multiphase, multicomponent mixture,
to denote the mole fraction of component
i in phase
j, and
to denote the mole fraction of component
i in the mixture.
is the total number of components and
is the total number of phases. The equilibrium ratio of component
i between phase
and phase
is defined as
. The multiphase Rachford-Rice equations that are derived from mass balances are
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
(
,
) and mole fractions
, Equation (1) is solved to give
through
. Then, the compositions of phases are obtained using
In this study, we introduce
to simply Equation (1) into a set of fractional equations
The values of vary between −1 and ∞.
Although Equation (4) appears simple, formulating a robust numerical solution is rather complex, because of the presence of
hyperplanes
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,
and
, defines a hypertetrahedron in the
space, and this hypertetrahedron cannot be dissected by any of the hyperplanes. This property ensures that
in Equation (4) are continuously differentiable in the hypervolume
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
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
. Our method began with the recognition that since
must be monotonically decreasing in the direction of increasing
. Because of this property, when
is increased to approach a singular hyperplane while all the other variables are kept constant,
. When
is decreased to approach a singular hyperplane while all the other variables are kept constant,
. Additionally, when
while all other variables are kept finite,
decreases and approaches zero; when
while all other variables are kept finite,
increases and approaches zero.
Figure 1 shows, qualitatively for a three-phase system, the monotonicity of
and
to
and
.
The monotonicity of to , 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, . Here, the subscripts to denote the number of iterative substitutions. We begin the first iteration by holding the values of through constant and seek to satisfy . Once is identified, we hold and through constant and seek to satisfy . This process continues till is solved and subscripts of all are updated to level “1”, which finishes the first iteration. The second iteration begins by holding the values of through constant and seek to satisfy , and ends when all the subscripts of are updated to level “2”. Iteration would then begin at level “3”, and continue till all are converged.
With this successive-substitution strategy, solution of Equation (4) is reduced to successive solutions of single-variable equation . Recognizing that is monotonic between its poles, we used the following method to solve :
Evaluate . Here, m refers to the level of iteration and in the square bracket is the variable that needs to be updated to level “m” such that .
Determine the direction along which to vary . If , should be increased. If, however, , should be decreased.
Check whether a solution exists along the direction of increasing (or decreasing) . This is achieved by performing a line search to see whether increasing (or decreasing) would lead to an intersection with a singular hyperplane. If needs to be increased and increasing generates intersections with singular hyperplanes, one solution must exist between the current and the nearest intersection, because, as approaches the intersection, approaches . If needs to be increased, but there is no intersection with singular hyperplanes along the direction of increasing , there is no solution to and the calculation stops. Similar criteria apply when needs to be decreased.
When solution exists, vary to find . Because there is one and only one solution between and the nearest intersection with singular hyperplanes, a bi-section method is used. Testing begins at the midpoint between and the nearest intersection. The interval that contains the solution is continuously halved until 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:
,
. 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
,
. 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
. Correlation in Hibbard and Schalla [
36] suggests that mole fraction of water in the oil phase can be set to
. The calculated compositions of the three phases and the equilibrium ratios are listed in
Table 1, together with their corresponding
.
This example was solved using a mixture composition of , and . The initial guess was . The case converged within five iterations with convergence criterion .
Figure 2 shows the path of convergence from the initial guess. Solution for
,
, and
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
and
. 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
and
. The dotted lines in
Figure 2 mark the locations where
and
, 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
, 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
. Assume that
is the result of the
-th iteration and
is the result of the
-th iteration, which is a successive-substitution. If the
-th iteration is a Newton step, we recommend not to use
to start the Newton step, but a combination of
and
. We specifically used
to start the Newton step after successive-substitution. Geometrically,
is the average of the points along the path of successive-substitution from
to
. As each of these points makes one and only one
zero,
calculated from Equation (8) is guaranteed to separate from surfaces
, 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 10
10. In all of the methods, the convergence criterion set on the difference between vectors
and
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 × 10
11 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 × 10
11. 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.