Rachford-Rice Equation Solvers (chemicals.rachford_rice)

This module contains functions for solving the Rachford-Rice Equation. This is used to solve ideal flashes, and is the inner loop of the sequential-substitution flash algorithm. It is not used by full newton-algorithms. The sequential-substitution is normally recommended because it does not suffer from the ~N^3 behavior of solving a matrix.

For reporting bugs, adding feature requests, or submitting pull requests, please use the GitHub issue tracker.

Two Phase - Interface

chemicals.rachford_rice.flash_inner_loop(zs, Ks, method=None, guess=None, check=False)

This function handles the solution of the inner loop of a flash calculation, solving for liquid and gas mole fractions and vapor fraction based on specified overall mole fractions and K values. As K values are weak functions of composition, this should be called repeatedly by an outer loop. Will automatically select an algorithm to use if no method is provided. Should always provide a solution.

The automatic algorithm selection will try an analytical solution, and use the Rachford-Rice method if there are 6 or more components in the mixture.

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

guessfloat, optional

Optional initial guess for vapor fraction, [-]

checkbool, optional

Whether or not to check the K values to ensure a positive-composition solution exists, [-]

Returns

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Other Parameters

methodstr, optional

The method name to use. Accepted methods are ‘Analytical’, ‘Rachford-Rice (Secant)’, ‘Rachford-Rice (Newton-Raphson)’, ‘Rachford-Rice (Halley)’, ‘Rachford-Rice (NumPy)’, ‘Leibovici and Nichita 2’, ‘Rachford-Rice (polynomial)’, and ‘Li-Johns-Ahmadi’. All valid values are also held in the list flash_inner_loop_methods.

Notes

A total of eight methods are available for this function. They are:

• ‘Analytical’, an exact solution derived with SymPy, applicable only only to mixtures of two, three, or four components

• ‘Rachford-Rice (Secant)’, ‘Rachford-Rice (Newton-Raphson)’, ‘Rachford-Rice (Halley)’, or ‘Rachford-Rice (NumPy)’, which numerically solves an objective function described in Rachford_Rice_solution.

• ‘Leibovici and Nichita 2’, a transformation of the RR equation described in Rachford_Rice_solution_LN2.

• ‘Li-Johns-Ahmadi’, which numerically solves an objective function described in Li_Johns_Ahmadi_solution.

• ‘Leibovici and Neoschil’, which numerically solves an objective function described in Rachford_Rice_solution_Leibovici_Neoschil.

Examples

>>> flash_inner_loop(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.6907302627738, [0.3394086969663, 0.3650560590371, 0.29553524399648], [0.571903654388, 0.27087159580558, 0.1572247498061])

chemicals.rachford_rice.flash_inner_loop_methods(N)

Return all methods able to solve the Rachford-Rice equation for the specified number of components.

Parameters

Nint

Number of components, [-]

Returns

methodslist[str]

Methods which can be used to solve the Rachford-rice equation

See also

flash_inner_loop

chemicals.rachford_rice.flash_inner_loop_all_methods = ('Analytical', 'Rachford-Rice (Secant)', 'Rachford-Rice (Newton-Raphson)', 'Rachford-Rice (Halley)', 'Rachford-Rice (NumPy)', 'Li-Johns-Ahmadi', 'Rachford-Rice (polynomial)', 'Leibovici and Nichita 2', 'Leibovici and Neoschil')

Tuple of method name keys. See the flash_inner_loop for the actual references

Two Phase - Implementations

chemicals.rachford_rice.Rachford_Rice_solution(zs, Ks, fprime=False, fprime2=False, guess=None)

Solves the objective function of the Rachford-Rice flash equation [1]. Uses the method proposed in [2] to obtain an initial guess.

$$\sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

fprimebool, optional

Whether or not to use the first derivative of the objective function in the solver (Newton-Raphson is used) or not (secant is used), [-]

fprime2bool, optional

Whether or not to use the second derivative of the objective function in the solver (parabolic Halley`s method is used if True) or not, [-]

guessfloat, optional

Optional initial guess for vapor fraction, [-]

Returns

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

The initial guess is the average of the following, as described in [2].

$$\left(\frac{V}{F}\right)_{min} = \frac{(K_{max}-K_{min})z_{of\;K_{max}} - (1-K_{min})}{(1-K_{min})(K_{max}-1)}$$

$$\left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}$$

Another algorithm for determining the range of the correct solution is given in [3]; [2] provides a narrower range however. For both cases, each guess should be limited to be between 0 and 1 as they are often negative or larger than 1.

$$\left(\frac{V}{F}\right)_{min} = \frac{1}{1-K_{max}}$$

$$\left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}$$

If the newton method does not converge, a bisection method (brenth) is used instead. However, it is somewhat slower, especially as newton will attempt 50 iterations before giving up.

In all benchmarks attempted, secant method provides better performance than Newton-Raphson or parabolic Halley`s method. This may not be generally true; but it is for Python and SciPy’s implementation. They are implemented for benchmarking purposes.

The first and second derivatives are:

$$\frac{d \text{ obj}}{d \frac{V}{F}} = \sum_i \frac{-z_i(K_i-1)^2} {(1 + \frac{V}{F}(K_i-1))^2}$$

$$\frac{d^2 \text{ obj}}{d (\frac{V}{F})^2} = \sum_i \frac{2z_i(K_i-1)^3} {(1 + \frac{V}{F}(K_i-1))^3}$$

References

1

Rachford, H. H. Jr, and J. D. Rice. “Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium.” Journal of Petroleum Technology 4, no. 10 (October 1, 1952): 19-3. doi:10.2118/952327-G.

2(1,2,3)

Li, Yinghui, Russell T. Johns, and Kaveh Ahmadi. “A Rapid and Robust Alternative to Rachford-Rice in Flash Calculations.” Fluid Phase Equilibria 316 (February 25, 2012): 85-97. doi:10.1016/j.fluid.2011.12.005.

3

Whitson, Curtis H., and Michael L. Michelsen. “The Negative Flash.” Fluid Phase Equilibria, Proceedings of the Fifth International Conference, 53 (December 1, 1989): 51-71. doi:10.1016/0378-3812(89)80072-X.

Examples

>>> Rachford_Rice_solution(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.6907302627738544, [0.33940869696634357, 0.3650560590371706, 0.2955352439964858], [0.5719036543882889, 0.27087159580558057, 0.15722474980613044])

chemicals.rachford_rice.Rachford_Rice_solution_LN2(zs, Ks, guess=None)

Solves the a objective function for the Rachford-Rice flash equation according to the Leibovici and Nichita (2010) transformation (method 2). This transformation makes the only zero of the function be the desired one. Consequently, higher-order methods may be used to solve this equation. Halley’s (second derivative) method is found to be the best; typically needing ~50% fewer iterations than the RR formulation with Secant method.

$$H(y) = \sum_i^n \frac{z_i}{\lambda - c_i} = 0$$

$$\lambda = c_k + \frac{c_{k+1} - c_k}{1 + e^{-y}}$$

$$c_i = \frac{1}{1-K_i}$$

$$c_{k} = \left(\frac{V}{F}\right)_{min}$$

$$c_{k+1} = \left(\frac{V}{F}\right)_{max}$$

Note the two different uses of c in the above equation, confusingly given in [1]. lambda is the vapor fraction.

Once the equation has been solved for y, the vapor fraction can be calculated outside the solver.

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

guessfloat, optional

Optional initial guess for vapor fraction, [-]

Returns

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

The initial guess is the average of the following, as described in [2].

$$\left(\frac{V}{F}\right)_{min} = \frac{(K_{max}-K_{min})z_{of\;K_{max}} - (1-K_{min})}{(1-K_{min})(K_{max}-1)}$$

$$\left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}$$

The first and second derivatives are derived with sympy as follows:

>>> from sympy import * 
>>> VF_min, VF_max, ai, ci, y = symbols('VF_min, VF_max, ai, ci, y') 
>>> V_over_F = (VF_min + (VF_max - VF_min)/(1 + exp(-y))) 
>>> F = ai/(V_over_F - ci) 
>>> terms = [F, diff(F, y), diff(F, y, 2)] 
>>> cse(terms, optimizations='basic') 

Some helpful information about this transformation can also be found in [3].

References

1

Leibovici, Claude F., and Dan Vladimir Nichita. “Iterative Solutions for ∑iaiλ-ci=1 Equations.” Chemical Engineering Research and Design 88, no. 5 (May 1, 2010): 602-5. https://doi.org/10.1016/j.cherd.2009.10.012.

2

Li, Yinghui, Russell T. Johns, and Kaveh Ahmadi. “A Rapid and Robust Alternative to Rachford-Rice in Flash Calculations.” Fluid Phase Equilibria 316 (February 25, 2012): 85-97. doi:10.1016/j.fluid.2011.12.005.

3

Billingsley, D. S. “Iterative Solution for ∑iaiλ-ci Equations.” Computers & Chemical Engineering 26, no. 3 (March 15, 2002): 457-60. https://doi.org/10.1016/S0098-1354(01)00767-0.

Examples

>>> Rachford_Rice_solution_LN2(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.6907302627738, [0.3394086969663, 0.3650560590371, 0.29553524399648], [0.571903654388, 0.27087159580558, 0.1572247498061])

chemicals.rachford_rice.Li_Johns_Ahmadi_solution(zs, Ks, guess=None)

Solves the objective function of the Li-Johns-Ahmadi flash equation. Uses the method proposed in [1] to obtain an initial guess.

$$0 = 1 + \left(\frac{K_{max}-K_{min}}{K_{min}-1}\right)x_{max} + \sum_{i=2}^{n-1} \frac{K_i-K_{min}}{K_{min}-1} \left[\frac{z_i(K_{max}-1)x_{max}}{(K_i-1)z_{max} + (K_{max}-K_i)x_{max}}\right]$$

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

Returns

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

The initial guess is the average of the following, as described in [1]. Each guess should be limited to be between 0 and 1 as they are often negative or larger than 1. max refers to the corresponding mole fractions for the species with the largest K value.

$$\left(\frac{1-K_{min}}{K_{max}-K_{min}}\right)z_{max}\le x_{max} \le \left(\frac{1-K_{min}}{K_{max}-K_{min}}\right)$$

If the newton method does not converge, a bisection method (brenth) is used instead. However, it is somewhat slower, especially as newton will attempt 50 iterations before giving up.

This method does not work for problems of only two components. K values are sorted internally. Has not been found to be quicker than the Rachford-Rice equation.

References

1(1,2)

Li, Yinghui, Russell T. Johns, and Kaveh Ahmadi. “A Rapid and Robust Alternative to Rachford-Rice in Flash Calculations.” Fluid Phase Equilibria 316 (February 25, 2012): 85-97. doi:10.1016/j.fluid.2011.12.005.

Examples

>>> Li_Johns_Ahmadi_solution(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.6907302627738544, [0.33940869696634357, 0.3650560590371706, 0.2955352439964858], [0.5719036543882889, 0.27087159580558057, 0.15722474980613044])

chemicals.rachford_rice.Rachford_Rice_solution_Leibovici_Neoschil(zs, Ks, guess=None)

Solves the objective function of the Rachford-Rice flash equation as modified by Leibovici and Neoschil. This modification helps convergence near the vapor fraction boundaries only; it slows convergence in other regions.

$$\left(\frac{V}{F} - \alpha_L\right)\left(\alpha_R - \frac{V}{F}\right) \sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

$$\alpha_L = - \frac{1}{K_{max} - 1}$$

$$\alpha_R = \frac{1}{1 - K_{min}}$$

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

guessfloat, optional

Optional initial guess for vapor fraction, [-]

Returns

L_over_Ffloat

Liquid fraction solution [-]

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

The initial guess is the average of the following.

$$\left(\frac{V}{F}\right)_{min} = \frac{(K_{max}-K_{min})z_{of\;K_{max}} - (1-K_{min})}{(1-K_{min})(K_{max}-1)}$$

$$\left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}$$

References

1

Leibovici, ClaudeF., and Jean Neoschil. “A New Look at the Rachford-Rice Equation.” Fluid Phase Equilibria 74 (July 15, 1992): 303-8. https://doi.org/10.1016/0378-3812(92)85069-K.

Examples

>>> Rachford_Rice_solution_Leibovici_Neoschil(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.3092697372261, 0.69073026277385, [0.339408696966343, 0.36505605903717, 0.29553524399648], [0.57190365438828, 0.270871595805580, 0.157224749806130])

chemicals.rachford_rice.Rachford_Rice_solution_polynomial(zs, Ks)

Solves the Rachford-Rice equation by transforming it into a polynomial, and then either analytically calculating the roots, or, using the known range the correct root is in, numerically solving for the correct polynomial root. The analytical solutions are used for N from 2 to 4.

Uses the method proposed in [2] to obtain an initial guess when solving the polynomial for the root numerically.

$$\sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

Warning

: Using this function with more than 20 components is likely to crash Python! This model does not work well with many components!

This method, developed first in [3] and expanded in [1], is clever but of little use for large numbers of components.

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

Returns

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

This approach has mostly been ignored by academia, despite some of its advantages.

The initial guess is the average of the following, as described in [2].

$$\left(\frac{V}{F}\right)_{min} = \frac{(K_{max}-K_{min})z_{of\;K_{max}} - (1-K_{min})}{(1-K_{min})(K_{max}-1)}$$

$$\left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}$$

If the newton method does not converge, a bisection method (brenth) is used instead. However, it is somewhat slower, especially as newton will attempt 50 iterations before giving up.

This method could be speed up somewhat for N <= 4; the checks for the vapor fraction range are not really needed.

References

1

Weigle, Brett D. “A Generalized Polynomial Form of the Objective Function in Flash Calculations.” Pennsylvania State University, 1992.

2(1,2)

Li, Yinghui, Russell T. Johns, and Kaveh Ahmadi. “A Rapid and Robust Alternative to Rachford-Rice in Flash Calculations.” Fluid Phase Equilibria 316 (February 25, 2012): 85-97. doi:10.1016/j.fluid.2011.12.005.

3

Warren, John H. “Explicit Determination of the Vapor Fraction in Flash Calculations.” Pennsylvania State University, 1991.

Examples

>>> Rachford_Rice_solution_polynomial(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.6907302627738543, [0.33940869696634357, 0.3650560590371706, 0.2955352439964858], [0.5719036543882889, 0.27087159580558057, 0.15722474980613044])

Two Phase - High-Precision Implementations

chemicals.rachford_rice.Rachford_Rice_solution_mpmath(zs, Ks, dps=200, tol=1e-100)

Solves the the Rachford-Rice flash equation using numerical root-finding to a high precision using the mpmath library.

$$\sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

dpsint, optional

Number of decimal places to use in the intermediate values of the calculation, [-]

tolfloat, optional

The tolerance of the solver used in mpmath, [-]

Returns

L_over_Ffloat

Liquid fraction solution [-]

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

This function is written solely for development purposes with the aim of returning bit-accurate solutions.

Note that the liquid fraction is also returned; it is insufficient to compute it as $\frac{L}{F} = 1 - \frac{V}{F}$.

Examples

>>> Rachford_Rice_solution_mpmath(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.3092697372261456, 0.6907302627738544, [0.33940869696634357, 0.3650560590371706, 0.29553524399648584], [0.5719036543882889, 0.27087159580558057, 0.15722474980613046])
>>> Rachford_Rice_solution_mpmath(zs=[0.999999999999, 1e-12], Ks=[2.0, 1e-12])
(1e-12, 0.999999999999, [0.49999999999975003, 0.50000000000025], [0.9999999999995001, 5.0000000000025e-13])

chemicals.rachford_rice.Rachford_Rice_solution_binary_dd(zs, Ks)

Solves the the Rachford-Rice flash equation for a binary system using double-double math. This increases the range in which the calculation can be performed accurately but does not totally eliminate error.

$$\sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

The analytical solution for a binary system is:

$$\frac{V}{F} = \frac{- K_{0} z_{0} - K_{1} z_{1} + z_{0} + z_{1}} {K_{0} K_{1} z_{0} + K_{0} K_{1} z_{1} - K_{0} z_{0} - K_{0} z_{1} - K_{1} z_{0} - K_{1} z_{1} + z_{0} + z_{1}}$$

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

Returns

L_over_Ffloat

Liquid fraction solution [-]

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Examples

This system with large volatility difference and a trace of a component shows a correct calculation. Try it out with other solvers for bad results!

>>> Rachford_Rice_solution_binary_dd(zs=[1E-27, 1.0], Ks=[1000000000000,0.1])
(1.000000000001, -1.0000000000009988e-12, [9.0000000000009e-13, 0.9999999999991], [0.90000000000009, 0.09999999999991001])

Note the limitations of this solver can be explored by comparing against Rachford_Rice_solution_mpmath. For example, with z0 of 1e-28 in the above example error creeps back in.

chemicals.rachford_rice.Rachford_Rice_solution_Leibovici_Neoschil_dd(zs, Ks, guess=None)

Solves the objective function of the Rachford-Rice flash equation as modified by Leibovici and Neoschil, using double-double precision math for maximum accuracy. For most cases, this function will return bit-for-bit accurate results; but there are pathological inputs where error still occurs.

$$\left(\frac{V}{F} - \alpha_L\right)\left(\alpha_R - \frac{V}{F}\right) \sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

$$\alpha_L = - \frac{1}{K_{max} - 1}$$

$$\alpha_R = \frac{1}{1 - K_{min}}$$

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

guessfloat, optional

Optional initial guess for vapor fraction, [-]

Returns

L_over_Ffloat

Liquid fraction solution [-]

V_over_Ffloat

Vapor fraction solution [-]

xslist[float]

Mole fractions of each species in the liquid phase, [-]

yslist[float]

Mole fractions of each species in the vapor phase, [-]

Notes

The initial guess is the average of the following.

$$\left(\frac{V}{F}\right)_{min} = \frac{(K_{max}-K_{min})z_{of\;K_{max}} - (1-K_{min})}{(1-K_{min})(K_{max}-1)}$$

$$\left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}$$

References

1

Leibovici, ClaudeF., and Jean Neoschil. “A New Look at the Rachford-Rice Equation.” Fluid Phase Equilibria 74 (July 15, 1992): 303-8. https://doi.org/10.1016/0378-3812(92)85069-K.

Examples

>>> Rachford_Rice_solution_Leibovici_Neoschil_dd(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
(0.3092697372261, 0.69073026277385, [0.339408696966343, 0.36505605903717, 0.29553524399648], [0.57190365438828, 0.270871595805580, 0.157224749806130])

Three Phase

chemicals.rachford_rice.Rachford_Rice_solution2(ns, Ks_y, Ks_z, beta_y=0.5, beta_z=1e-06)

Solves the two objective functions of the Rachford-Rice flash equation for a three-phase system. Initial guesses are required for both phase fractions, beta_y and beta_z. The Newton method is used, with an analytical Jacobian.

$$F_0 = \sum_i \frac{z_i (K_y -1)}{1 + \beta_y(K_y-1) + \beta_z(K_z-1)} = 0$$

$$F_1 = \sum_i \frac{z_i (K_z -1)}{1 + \beta_y(K_y-1) + \beta_z(K_z-1)} = 0$$

Parameters

nslist[float]

Overall mole fractions of all species (would be zs except that is conventially used for one of the three phases), [-]

Ks_ylist[float]

Equilibrium K-values of y phase to x phase, [-]

Ks_zlist[float]

Equilibrium K-values of z phase to x phase, [-]

beta_yfloat, optional

Initial guess for y phase (between 0 and 1), [-]

beta_zfloat, optional

Initial guess for z phase (between 0 and 1), [-]

Returns

beta_yfloat

Phase fraction of y phase, [-]

beta_zfloat

Phase fraction of z phase, [-]

xslist[float]

Mole fractions of each species in the x phase, [-]

yslist[float]

Mole fractions of each species in the y phase, [-]

zslist[float]

Mole fractions of each species in the z phase, [-]

Notes

The elements of the Jacobian are calculated as follows:

$$\frac{\partial F_0}{\partial \beta_y} = \sum_i \frac{-z_i (K_y -1)^2} {\left(1 + \beta_y(K_y-1) + \beta_z(K_z-1)\right)^2}$$

$$\frac{\partial F_1}{\partial \beta_z} = \sum_i \frac{-z_i (K_z -1)^2} {\left(1 + \beta_y(K_y-1) + \beta_z(K_z-1)\right)^2}$$

$$\frac{\partial F_1}{\partial \beta_y} = \sum_i \frac{\partial F_0} {\partial \beta_z} = \frac{-z_i (K_z -1)(K_y - 1)}{\left(1 + \beta_y(K_y-1) + \beta_z(K_z-1)\right)^2}$$

In general, the solution which Newton’s method converges to may not be the desired one, so further constraints are required.

Okuno’s method in [1] provides a polygonal region where the correct answer lies. It has not been implemented.

The Leibovici and Neoschil method [4] provides a method to compute/update the damping parameter, which is suposed to ensure convergence. It claims to be able to calculate the maximum damping factor for Newton’s method, if it tries to go out of bounds.

A custom region which is believed to be the same as that of Okuno is implemented instead - the region which ensures positive compositions for all compounds in all phases, but does not restrict the phase fractions to be between 0 and 1 or even positive.

With the convergence restraint, it is believed if a solution lies within (0, 1) for both variables, the correct solution will be converged to so long as the initial guesses are within the correct region.

Some helpful information has also been found in [2] and [3].

References

1

Okuno, Ryosuke, Russell Johns, and Kamy Sepehrnoori. “A New Algorithm for Rachford-Rice for Multiphase Compositional Simulation.” SPE Journal 15, no. 02 (June 1, 2010): 313-25. https://doi.org/10.2118/117752-PA.

2

Li, Zhidong, and Abbas Firoozabadi. “Initialization of Phase Fractions in Rachford-Rice Equations for Robust and Efficient Three-Phase Split Calculation.” Fluid Phase Equilibria 332 (October 25, 2012): 21-27. https://doi.org/10.1016/j.fluid.2012.06.021.

3

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

4

Leibovici, Claude F., and Jean Neoschil. “A Solution of Rachford-Rice Equations for Multiphase Systems.” Fluid Phase Equilibria 112, no. 2 (December 1, 1995): 217-21. https://doi.org/10.1016/0378-3812(95)02797-I.

Examples

>>> ns = [0.204322076984, 0.070970999150, 0.267194323384, 0.296291964579, 0.067046080882, 0.062489248292, 0.031685306730]
>>> Ks_y = [1.23466988745, 0.89727701141, 2.29525708098, 1.58954899888, 0.23349348597, 0.02038108640, 1.40715641002]
>>> Ks_z = [1.52713341421, 0.02456487977, 1.46348240453, 1.16090546194, 0.24166289908, 0.14815282572, 14.3128010831]
>>> Rachford_Rice_solution2(ns, Ks_y, Ks_z, beta_y=.1, beta_z=.6)
(0.6868328915094766, 0.06019424397668606, [0.1712804659711611, 0.08150738616425436, 0.1393433949193188, 0.20945175387703213, 0.15668977784027893, 0.22650123851718007, 0.015225982711774586], [0.21147483364299702, 0.07313470386530294, 0.31982891387635903, 0.33293382568889657, 0.036586042443791586, 0.004616341311925655, 0.02142533917172731], [0.26156812278601893, 0.00200221914149187, 0.20392660665189805, 0.2431536850887592, 0.03786610596908295, 0.03355679851539993, 0.21792646184834918])

N Phase

chemicals.rachford_rice.Rachford_Rice_solutionN(ns, Ks, betas)

Solves the (phases -1) objectives functions of the Rachford-Rice flash equation for an N-phase system. Initial guesses are required for all phase fractions except the last. The Newton method is used, with an analytical Jacobian.

Parameters

nslist[float]

Overall mole fractions of all species, [-]

Kslist[list[float]]

Equilibrium K-values of all phases with respect to the x (reference) phase, [-]

betaslist[float]

Phase fraction initial guesses only for the first N - 1 phases; each value corresponds to the phase fraction of each set of the K values; if a phase fraction is specified for the last phase as well, it is ignored [-]

Returns

betaslist[float]

Phase fractions of all of the phases; one each for each K value set given, plus the reference phase phase fraction [-]

compositionslist[list[float]]

Mole fractions of each species in each phase; order each phase in the same order as the K values were provided, and then the x phase last, which was the reference phase [-]

Notes

This algorithm has been used without issue for 4 and 5 phase flashes.

Some helpful information was found in [1], although this method does not follow it exactly.

References

1

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

Examples

>>> ns = [0.204322076984, 0.070970999150, 0.267194323384, 0.296291964579, 0.067046080882, 0.062489248292, 0.031685306730]
>>> Ks_y = [1.23466988745, 0.89727701141, 2.29525708098, 1.58954899888, 0.23349348597, 0.02038108640, 1.40715641002]
>>> Ks_z = [1.52713341421, 0.02456487977, 1.46348240453, 1.16090546194, 0.24166289908, 0.14815282572, 14.3128010831]
>>> Rachford_Rice_solutionN(ns, [Ks_y, Ks_z], [.1, .6])
([0.6868328915094767, 0.06019424397668605, 0.25297286451383727], [[0.21147483364299702, 0.07313470386530294, 0.3198289138763589, 0.33293382568889657, 0.03658604244379159, 0.004616341311925657, 0.02142533917172731], [0.26156812278601893, 0.00200221914149187, 0.203926606651898, 0.2431536850887592, 0.03786610596908296, 0.033556798515399944, 0.21792646184834918], [0.1712804659711611, 0.08150738616425436, 0.13934339491931877, 0.20945175387703213, 0.15668977784027896, 0.22650123851718015, 0.015225982711774586]])

Two Phase Utility Functions

chemicals.rachford_rice.Rachford_Rice_polynomial(zs, Ks)

Transforms the Rachford-Rice equation into a polynomial and returns its coefficients. A spelled-out solution is used for N from 2 to 5, derived with SymPy and optimized with the common sub expression approach.

Warning

For large numbers of components (>20) this model performs terribly, though with future optimization it may be possible to have better performance.

$$\sum_{i=1}^N z_i C_i\left[ \Pi_{j\ne i}^N \left(1 + \frac{V}{F} C_j\right)\right] = 0$$

$$C_i = K_i - 1.0$$

Once the above calculation is performed, it must be rearranged into polynomial form.

Parameters

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

Returns

coeffsfloat

Coefficients, with earlier coefficients corresponding to higher powers, [-]

Notes

Explicit calculations for any degree can be obtained with SymPy, changing N as desired:

>>> from sympy import * 
>>> N = 4
>>> Cs = symbols('C0:' + str(N)) 
>>> zs = symbols('z0:' + str(N)) 
>>> alpha = symbols('alpha') 
>>> tot = 0
>>> for i in range(N): 
...     mult_sum = 1
>>> for j in range(N): 
...     if j != i:
...         mult_sum *= (1 + alpha*Cs[j])
...     tot += zs[i]*Cs[i]*mult_sum

poly_expr = poly(expand(tot), alpha) coeff_list = poly_expr.all_coeffs() cse(coeff_list, optimizations=’basic’)

[1] suggests a matrix-math based approach for solving the model, but that has not been performed here. [1] also has explicit equations for up to N = 7 to derive the coefficients.

The general form was derived to be slightly different than that in [1], but is confirmed to also be correct as it matches other methods for solving the Rachford-Rice equation. [2] has similar information to [1].

The first coefficient is always 1.

The approach is also discussed in [3], with one example.

References

1(1,2,3,4)

Weigle, Brett D. “A Generalized Polynomial Form of the Objective Function in Flash Calculations.” Pennsylvania State University, 1992.

2

Warren, John H. “Explicit Determination of the Vapor Fraction in Flash Calculations.” Pennsylvania State University, 1991.

3

Monroy-Loperena, Rosendo, and Felipe D. Vargas-Villamil. “On the Determination of the Polynomial Defining of Vapor-Liquid Split of Multicomponent Mixtures.” Chemical Engineering Science 56, no. 20 (October 1, 2001): 5865-68. https://doi.org/10.1016/S0009-2509(01)00267-6.

Examples

>>> Rachford_Rice_polynomial(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
[1.0, -3.6926529966760824, 2.073518878815093]

chemicals.rachford_rice.Rachford_Rice_flash_error(V_over_F, zs, Ks)

Calculates the objective function of the Rachford-Rice flash equation. This function should be called by a solver seeking a solution to a flash calculation. The unknown variable is V_over_F, for which a solution must be between 0 and 1.

$$\sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

Parameters

V_over_Ffloat

Vapor fraction guess [-]

zslist[float]

Overall mole fractions of all species, [-]

Kslist[float]

Equilibrium K-values, [-]

Returns

errorfloat

Deviation between the objective function at the correct V_over_F and the attempted V_over_F, [-]

Notes

The derivation is as follows:

$$F z_i = L x_i + V y_i$$

$$x_i = \frac{z_i}{1 + \frac{V}{F}(K_i-1)}$$

$$\sum_i y_i = \sum_i K_i x_i = 1$$

$$\sum_i(y_i - x_i)=0$$

$$\sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0$$

This objective function was proposed in [1].

References

1

Rachford, H. H. Jr, and J. D. Rice. “Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium.” Journal of Petroleum Technology 4, no. 10 (October 1, 1952): 19-3. doi:10.2118/952327-G.

Examples

>>> Rachford_Rice_flash_error(0.5, zs=[0.5, 0.3, 0.2],
... Ks=[1.685, 0.742, 0.532])
0.04406445591174976

Numerical Notes

For the two-phase problem, there are the following ways of computing the vapor and liquid mole fractions once the vapor fraction and liquid fraction has been computed:

The most commonly shown expression is:

$$x_i = \frac{z_i}{1 + \frac{V}{F}(K_i-1)}$$

This can cause numerical issues when $K_i$ is near 1. It also shows issues near $\frac{V}{F}(K_i-1) = -1$.

Another expression which avoids the second issue is

$$x_i = \frac{z_i}{\frac{L}{F} + (1 - \frac{L}{F})K_i}$$

Much like the other expression above this numerical issues but at different conditions: $\frac{L}{F} = 1$ and $\frac{L}{F} = -(1 - \frac{L}{F})K_i$.

One more expression using both liquid and vapor fraction is:

$$x_i = \frac{z_i}{K_i\frac{V}{F} + \frac{L}{F} }$$

This expression only has one problematic area: $K_i\frac{V}{F} = \frac{L}{F}$. Preferably, this is computed with a fused-multiply-add operation.

Another expression which flips the K value into the liquid form and swaps the vapor fraction for the liquid fraction in-line is as follows

$$x_i = \frac{\frac{z_i}{K_i}} {\frac{\frac{L}{F}}{K_i} + \frac{V}{F}}$$

This also has numerical problems when $-\frac{\frac{L}{F}}{K_i} = \frac{V}{F}$.

Even when computing a solution with high precision such as with mpmath, the resulting compositions and phase fractions may fail basic tests. In the following case, a nasty problem has a low-composition but relatively volatile last component. Mathematically, $1 = \frac{\frac{L}{F} x_i + \frac{V}{F} y_i}{z_i}$. This is true for all components except the last one in this case, where significant error exists.

>>> zs = [0.004632150100959984, 0.019748784459594933, 0.0037494212674659875, 0.0050492815033649835, 7.049818284201636e-05, 0.019252941309184937, 0.022923068733233923, 0.02751809363371991, 0.044055273670258854, 0.026348159124199914, 0.029384949788372902, 0.022368938441593926, 0.03876345111451487, 0.03440715821883388, 0.04220510198067186, 0.04109191458414686, 0.031180945124537895, 0.024703227642798916, 0.010618543295340965, 0.043262442161003854, 0.006774922650311977, 0.02418090788262392, 0.033168278052077886, 0.03325881573680989, 0.027794535589044905, 0.00302091746847699, 0.013693571363003955, 0.043274465132840854, 0.02431371852108292, 0.004119055065872986, 0.03314056562191489, 0.03926511182895087, 0.0305068048046159, 0.014495317922126952, 0.03603737707409988, 0.04346278949361786, 0.019715052322446934, 0.028565255195219907, 0.023343683279902924, 0.026532427286078915, 2.0833722372767433e-06]
>>> Ks = [0.000312001984979, 0.478348350355814, 0.057460349529956, 0.142866526725442, 0.186076915390803, 1.67832923245552, 0.010784509466239, 0.037204384948088, 0.005359146955631, 2.41896552551221, 0.020514598049597, 0.104545054017411, 2.37825397780443, 0.176463709057649, 0.000474240879865, 0.004738042026669, 0.02556030236928, 0.00300089652604, 0.010614774675069, 1.75142303167203, 1.47213647779132, 0.035773024794854, 4.15016401471676, 0.024475125100923, 0.00206952065986, 2.09173484409107, 0.06290795470216, 0.001537212006245, 1.16935817509767, 0.001830422812888, 0.058398776367331, 0.516860928072656, 1.03039372722559, 0.460775800103578, 0.10980302936483, 0.009883724220094, 0.021938589630783, 0.983011657214417, 0.01978995396409, 0.204144939961852, 14.0521979447538]
>>> LF, VF, xs, ys = Rachford_Rice_solution_mpmath(zs=zs, Ks=Ks)
>>> (LF*xs[-1] + VF*ys[-1])/zs[-1]
1.0000000000028162