Module Goal: To use cubic equations of state for the description of phase equilibrium.
Module Objective: To integrate the knowledge of qualitative phase behavior and the framework for quantification to actually build a phase equilibrium model.
Gas and liquid co-existence is common to all petroleum and natural gas applications, as we have mentioned already. Additionally, in our discussions presented in Modules 12 and 13, we tried to mathematically model the problem of vapor-liquid co-existence in equilibrium. As we recall, at that time we concluded that we were lacking tools necessary to complete the definition of the problem. Then we went on and concentrated on building those thermodynamic tools which we needed. Now we are ready to couple the concepts reviewed in Modules 14 through 16 (“Thermodynamic Tools”) with those reviewed in Modules 12 and 13. In order for us to do that, we will need to apply our knowledge of equations of state (studied in Modules 6 through 11).
Once the mathematics of vapor-liquid equilibrium is understood, we will be able to solve for a variety of applications. They will be discussed in the last modules of this course. One of those applications will be to generate any sort of phase diagrams, like the ones that we studied in our earlier discussions.
As you see, everything that we have studied in this course is completely interlinked. That’s one of the main ideas that I want you to glean from these discussions.
Consider a liquid-vapor in equilibrium. As we have discussed previously, a condition for equilibrium is that the chemical potential of each component in both phases are equal, thus:
We showed that this is equivalent to:
This is, for a system to be in equilibrium, the fugacity of each component in each of the phases must be equal as well. The fugacity of a component in a mixture can be expressed in terms of the fugacity coefficient. Therefore, the fugacity of a component in either phase can be written as:
Introducing (17.3) into (17.2),
This equilibrium condition can be written in terms of the equilibrium ratio ${K}_{i}={y}_{i}l{x}_{i}$ , to get:
Do you recall the problem at the end of Module 13? At that point we needed a more reliable way to calculate the equilibrium ratios that showed up in the Rachford-Rice objective function. We demonstrated that once we know all values of K_{i}’s, the problem of vapor-liquid equilibrium is reduced to solving the Rachford-Rice objective function, using the Newton-Raphson Procedure.
We can now calculate equilibrium ratios, using (17.5), in terms of fugacity coefficients. We also know that we have an analytic expression for the calculation of fugacity coefficients via EOS — this was shown in the last section of the previous module. This is why we call this module “Vapor Liquid Equilibrium via EOS.”
Is this the end to our problems? Not quite. Take a look at the expression for fugacity coefficients in mixtures both for SRK EOS and PR EOS. It is clear that they are functions of the pressure, temperature, and composition of the phases:
Do we know the composition of the phases “x_{i}”, “y_{i}” in advance? In a typical flash problem, we are given pressure, temperature and overall composition (z_{i}). What do we want to know? How much gas, how much liquid, and the compositions of the phases: ${\alpha}_{g},\text{}{\alpha}_{l},{\text{y}}_{i},{\text{x}}_{i}$ . So, we do not know those compositions in advance; therefore, as it stands, we cannot calculate (17.6) or (17.5). Thus far, it seems that the flash problem is unsolvable.
If we are bold enough, we could try to overcome this problem by “guessing” those compositions, and proceed by solving (17.6) and (17.5). With this “rough” estimate for Ki’s, we could solve for “${\alpha}_{g}$ ” with the procedure outlined in Module 13 (“Objective Function and Newton-Raphson Procedure”). Once “${\alpha}_{g}$” is known, we could back calculate the compositions of the phases using equations (12.7) and (12.11). If we were correct, those compositions would match each other (the “guessed” ones with respect to the “back-calculated”). More than likely, this would not happen, and we would have to make a new “guess.” This is, fundamentally, an iterative procedure. Although this is not what we do, it does illustrate that this problem is solvable by implementing the appropriate iterative scheme.
In equations (17.4) and (17.5), the fugacity of the liquid and vapor phases were computed in terms of the fugacity coefficient. Hence, this method of expressing the equilibrium criteria is known as the dual-fugacity coefficient method. For the sake of completeness, it is necessary to indicate that the fugacity of a component in a mixture can also be expressed in terms of a thermodynamic concept called the activity coefficient. While the fugacity coefficient is seen as a measure of the deviation of behavior with respect to the ideal gas model, the activity coefficient measures the deviation of behavior with respect to the ideal liquid model. This approach is called the dual-activity coefficient method, in which both liquid and vapor phase fugacities are expressed in terms of the activity coefficient and substituted into the equilibrium criteria in (17.2). A mixed activity coefficient-fugacity coefficient method can be also devised by expressing the liquid phase fugacities in terms of activity coefficients and the vapor phase fugacities in terms of fugacity coefficients. Each of the aforementioned methods for the calculation of phase equilibria has its advantages and disadvantages. The dual-fugacity-coefficient method is simpler both conceptually and computationally, but if the equation of state does not predict liquid and vapor densities well, the results may be inaccurate. The activity coefficient method can be more accurate, but it is more complicated to implement. For the rest of the discussion, the dual-fugacity coefficient approach will be used.
We have essentially reduced the typical VLE problem to that of solving a system of non-linear algebraic equations. The objective function is the nucleus of the VLE calculation. Rather than solve the material balance equations in their elemental form, the use of an objective function will allow for a more stable iteration scheme. We have seen that three objective functions are available for this calculation:
Expression (17.7c) is the Rachford-Rice Objective Function; we have proven that its numerical solution is cleaner than that of equations (17.7a) and (17.7b). We also discussed that the key problem in solving for “${\alpha}_{g}$” is the need for values for the equilibrium ratios (K_{i}’s). With the aid of equations (17.5) and (17.6) these problems can be solved simultaneously by iteration techniques.
In the problem we have posed, we have (3n+1) unknowns:
And, ${\alpha}_{g}$
And we also have (3n+1) equations:
where P, T, and z_{i} are the values that we know. See how this is only a restatement of the original set of equations that contain the equilibrium conditions, material balances, and molar fraction constraints:
Since we have the same number of equations as unknowns, this system is determined and has a unique solution. To solve this system of simultaneous non-linear equations, there are two categories of solution techniques:
Newton-type methods for more than one unknown require finding (n x n) elements of a Jacobian matrix at every iteration step. This could be computationally expensive. Additionally, the Newton Raphson method requires a very good initial guess (i.e., close to the actual solution) to guarantee convergence to the true values. This is not always possible, especially at the start of the procedure.
The most popular method, and the easiest to implement, is the substitution type. However, substitution type methods can be quite slow for some conditions of interest. In these sort of cases, we either switch to a Newton-Raphson solver (that has a much better rate of convergence) or implement an acceleration procedure to the substitution method itself.
In a substitution-type method, we start with initial guesses for all of the unknowns and loop around the equations to obtain “better” approximations for each of them. We test the goodness of the solution at every time step by comparing the new, better approximation to the previous guess. If the correction is small under certain convergence criteria, the procedure is stopped and we use the results from the last iteration as the final answer.
As we discussed in the previous section, reliable values for the equilibrium constants (K_{i}’s) must be obtained before we can solve the Rachford-Rice Objective Function. Generally, first estimates for equilibrium constants are calculated by using Wilson’s empirical equation (equation 13.15). However, Wilson’s correlation yields only approximate values for equilibrium ratios. Here, we apply thermodynamic equilibrium considerations in order to obtain more reliable predictions for K_{i}’s.
Let us recall that for a system to be in equilibrium, any net transfer (of heat, momentum, or mass) must be zero. For this to be true, all potentials (temperature, pressure, and chemical potential) must be the same in all of the phases. Provided that the temperature and pressure of both phases are the same, a zero net transfer for all components in the mixture results when all chemical potentials are the same. A restatement of this: all fugacities for all components in each phase are equal. Since fugacity is a measure of the potential for transfer of a component between two phases, equal fugacities of a component in both phases results in zero net transfer. SSM (Successive Substitution Method) takes advantage of this fact. Equation (17.5) showed that the fugacity criterion for equilibrium allows writing the equilibrium ratios K_{i} as a function of fugacity coefficients as follows:
Equation (17.10) presupposes equality of fugacity (${f}_{i}^{V}={f}_{i}^{L}$), but during an iterative procedure, fugacities may not be equal (${f}_{i}^{V}\ne {f}_{i}^{L}$) until convergence conditions have been attained. Therefore, if fugacities are not equal (convergence has not been achieved), equation (17.10) becomes:
Using the above equation, a correction step can be formulated to improve the current values of K_{i}. The SSM-updating step is written as:
SSM updates all previous equilibrium ratios (K_{i}) using the fugacities predicted by the equation of state. This iteration method requires an initial estimation of K_{i}-values — for which Wilson’s correlation is used. It can easily be concluded that the convergence criteria will be satisfied whenever the fugacity ratios of all the components in the system are close to unity. Such condition is achieved when the following inequality is satisfied:
When the system is close to the critical point and fugacities are strongly composition-dependent, a slowing-down of the convergence rate of the SSM (Successive Substitution Method) is to be expected. In an attempt to avoid slow convergence problems, some methods have been proposed. Among the most popular are the Minimum Variable Newton Raphson (MVNR) Method and the Accelerated and Stabilized Successive Substitution Method (ASSM).
The ASSM is basically an accelerated version of the SSM procedure, and thus follows a similar theory. Such procedure is implemented to accelerate the calculation of K_{i}-values, especially in the region close to critical point where the use of the SSM alone will not be efficient. The ASSM technique was presented by Rinses et al. (1981) and consists of the following steps:
These criteria show that you have sufficient proximity to the conditions to ensure the efficiency of the method. Rr_{i} is the ratio of liquid fugacity to gas fugacity of the i-th component and ‘${\alpha}_{g}$’ is molar gas fraction of the two-phase system.
If the system satisfies ALL above criteria, the iteration technique is then switched from the SSM to the ASSM. Otherwise, SSM is used for the update of the K_{i}-values. The following expressions are used to update K_{i}-values in ASSM:
where ${\lambda}_{i}=\left[\left(R{r}_{i}^{old}-1\right)/\left(R{r}_{i}^{old}-R{r}_{i}^{new}\right)\right]$
In some cases, using a constant acceleration value of ${\lambda}_{i}=2$ is good enough.
Even though we are outlining Risnes et al.’s version of the accelerated Successive Substitution Method, there are several other published algorithms whose main purpose have also been to accelerate the successive substitution method. Fussel and Yanosik (1978), Michelsen (1982), and Mehra et al. (1983) are examples of such attempts. Risnes et al. version is the easiest and most straightforward to implement, but it is subjected to limitations.
Interestingly enough, one of the most difficult aspects of making VLE calculations may not be the two-phase splitting calculation itself, but knowing whether or not a mixture will actually split into two (or even more) phases for a pressure and temperature condition.
A single-phase detection routine has to be simultaneously introduced at this stage to detect whether the system is in a true single-phase condition at the given pressure and temperature or whether it will actually split into two-phases. Several approaches may be used here: the Bring-Back technique outlined by Risnes et al. (1981), and Phase Stability Criteria introduced by Michelsen (1982), among others. Here we describe Michelsen’s stability test.
Michelsen (1982) suggested creating a second-phase inside any given mixture to verify whether such a system is stable or not. It is the same idea behind the Bring-Back procedure (Risnes et al., 1981), but this test additionally provides straightforward interpretation for the cases where trivial solutions are found
(K_{i}’s —> 1). The test must be performed in two parts, considering two possibilities: the second phase can be either vapor-like or liquid-like. The outline of the method is described below, following the approach presented by Whitson and Brule (2000).
If a trivial solution is approached, stop the procedure.
If convergence has not been attained, use the new K-values and go back to step (b).
Create a liquid-like second phase,
Follow the previous steps by replacing equations (17.15), (17.16), (17.17), and (17.18) by (17.22), (17.23), (17.24), and (17.25) respectively.
The interpretation of the results of this method follows:
Answer the following problem, and submit your answer to the drop box in Canvas that has been created for this module.
Please note: