Rapid multi-component phase-split calculations using volume functions and reduction methods

We present a new family of fast and robust methods for the calculation of the vapor-liquid equilibrium at isobaric-isothermal (PT-flash), isochoric-isothermal (VT-flash), isenthalpic-isobaric (HP-flash), and isoenergetic-isochoric (UV-flash) conditions. The framework is provided by formulating phase-equilibrium conditions for multi-component mixtures in an effectively reduced space based on the molar specific value of the recently introduced volume function derived from the Helmholtz free energy. The proposed algorithmic implementation can fully exploit the optimum quadratic convergence of a Newton method with the analytical Jacobian matrix. This paper provides all required exact analytic expressions for the general cubic equation of state. Computational results demonstrate the effectivity and efficiency of the new methods. Compared to conventional methods, the proposed reduced-space iteration leads to a considerable speed-up as well as to improved robustness and better convergence behavior near the spinodal and coexistence curves of multi-component mixtures, where the preconditioning by the reduction method is most effective.


Introduction
Robust, computationally efficient and accurate phase splitting or flash calculations play a crucial role in many engineering disciplines, such as chemical-process and reservoir simulations. In Com-putational Fluid Dynamics (CFD) simulations of realistic multi-component vapor-liquid fluid flows, millions of phase equilibrium calculations are required every time step in the form of either the VT-flash or UV-flash, depending on the chosen formulation of the governing equations: The VTflash is needed in cases where the overall specific volume, temperature and composition are known, such as for the carbon dioxide injection into subsurface reservoirs 1,2 . Methods that solve the compressible Navier-Stokes equations based on the conservation laws for mass, linear momentum and total energy, such as applied for the simulation of the trans-critical vaporization of liquid fuels 3-6 , require a UV-flash, where the input is the overall specific internal energy, volume and composition.
The calculation of thermodynamic equilibrium properties of multi-component multi-phase mixtures typically consumes more than three quarters of the total computational time 7,8 and thus imposes severe limitation on the tractable space-time resolution or even the computational feasibility of such numerical simulations. At the same time, flash algorithms for CFD applications have to be fault tolerant and robust, because even a method that fails to converge only once in a billion will eventually spoil the entire simulation.
The simplest case and workhorse of most phase-equilibrium calculations is the so-called PT-flash, where the equilibrium pressure and temperature of the mixture are already given. Most methods for calculating the isobaric-isothermal equilibrium volume fractions and compositions follow the approach proposed by Michelsen 9,10 . For solving flash problems at conditions other than constant pressure and temperature, Michelsen 11 introduced an indirect method based on nesting simple and robust PT-flash calculations. For VT-flashes, for example, Michelsen's method aims to find the pressure at which the PT-flash results in the given total specific volume. This results in an optimization algorithm, in which the pressure is adjusted in the outer loop and a PT-flash is solved in the inner loop. Accordingly, UV-flashes are solved by a bi-variant optimization of temperature and pressure corresponding to the given internal energy and volume, which define the thermodynamic state, for example, in mass and energy conservative Navier-Stokes solvers 3 .
Nested algorithms based on the PT-flash are also attractive for mixtures with many components, because they offer the possibility of adopting reduction methods 12 , which provide a considerable speedup and, in addition, improve the robustness of the algorithm 13 . The first reduction method was introduced by Michelsen 12 , who found that the phase-splitting problem is fully defined by only three reduced parameters regardless of the number of components when all Binary Interaction Coefficients (BICs) are zero. Hendriks and Van Bergen 14 successfully generalized the method for cases with some non-zero BICs through an eigenvalue analysis of the binary interaction matrix.
Nichita and Graciaa 15 found a new set of reduced parameters for PT-flash calculations, for which they demonstrated a notable decrease in the number of iterations relative to previous reduction methods specifically near the phase boundary and the critical point.
Employing a direct VT-flash, on the other hand, could considerably reduce the computational time by eliminating the outer pressure iteration loop, provided that the method itself would be fast and robust enough. To this end, Mikyška and Firoozabadi 16 introduced an alternative formulation of the VT-flash problem based on a new thermodynamic function, the so-called volume function.
They solved the problem directly by a successive-substitution iteration (SSI) algorithm with nearly the same number of iterations as a conventional PT-flash based solver requires for one inner iteration loop. Recently, Jindrová and Mikyška 17 and Nichita 18 presented methods for solving the VT-flash problem via direct minimization of the total Helmholtz free energy. Cismondi et al. 19 directly included the pressure equality and volume constraint in a new algorithm very similar to the PTflash, and showed about 20% reduction in the computational time compared to Michelsen's nested optimization technique. However, for working fluids with a large number of components, these methods lead to a significantly stronger increase of the computational time compared to the nested approach that benefits from the quadratic Newton-Raphson convergence rate in the reduced space.
Extending the work of Mikyška and Firoozabadi 16 and Nichita and Graciaa 15 , we present a very fast and robust method for direct vapor-liquid phase-split calculations based on formulating phase equilibrium conditions in terms of the molar specific value of Mikyška and Firoozabadi's volume function (instead of fugacity coefficients) and a corresponding reduction method. This new formulation allows us to solve isothermal flashes (both PT and VT) directly and with the exact analytical Jacobian matrix, which results in optimum quadratic convergence of the Newton-Raphson method. Non-isothermal cases, such as UV and HP flashes, are solved through nested univariate optimization with the corresponding isothermal flash (PT for HP and VT for UV) and the readily available specific heat capacity at constant pressure (for HP-flash) or at constant volume (for UV-flash) as exact Jacobian. This paper is structured as follows: First, the mathematical description of the equilibrium problem is reformulated based on the molar specific values of the volume function for the vapor and liquid phases. Next, the classical reduction method is presented along with the derivation of new reduced parameters in the context of the new formulations, and all other thermodynamic relations required for non-isothermal flashes are also derived from the reduced parameters. Then, algorithms based on the Newton-Raphson method with the analytical Jacobian matrix for the direct solution of isothermal flashes and for the indirect solution of non-isothermal ones, are presented. Last but not least, the reliability and efficiency of proposed algorithms and its significantly improved computational performance compared to a recently published implementation for high-fidelity CFD simulations 3 will be demonstrated and discussed for different multi-component mixtures at various thermodynamic conditions.

Thermodynamic equilibrium formulation
According to the Gibbsian thermodynamics 20 , a multi-component system consisting of vapor and liquid phases is in equilibrium when the temperatures, pressures, and chemical potentials of phases are equal, that is, where T , p, and µ i are temperature, pressure, and chemical potential of component i = {1 . . . n} in a mixture with n components, and superscripts L and V refer to values of the liquid and vapor phases.
The pressures can be computed as a function of temperature, molar specific volume and composition of each phase using the general form of the cubic equation of states (EoS) where δ 1 and δ 2 are the two EoS parameters (see below), R is the universal gas constant, and v is the molar specific volume of the mixture. The energy and co-volume parameters a and b are usually computed using the van der Waals mixing rules in which z i is the mole fraction of the component i, and κ ij is the binary interaction coefficient between component i and j in the mixture.â i andb i are energy and co-volume parameters of the pure component i, which are obtained througĥ where T ci , p ci are critical temperature and pressure of the component i. The two constants Ω a and Ω b as well as the form of the function of c(ω i ), in which ω i is the acentric factor, depend of the selected cubic EoS: for instance, in the Peng-Robinson (PR) EoS: which result in Ω a = 0.45724, Ω b = 0.0778, and the functional c(ω i ) is: In the Soave-Redlich-Kwong (SRK) EoS, with δ 1 = 0 and δ 2 = 1, they are Ω a = 0.42748, Ω b = 0.08664, The equality of chemical potentials is typically expressed in terms of the K-factor (also named K-value or equilibrium ratio, which is the ratio of the mole fractions in the vapor (y) and liquid (x) phases), and the fugacity coefficient derived from the Gibbs free energy. The logarithmic form of this relation for a two-phase vapor-liquid mixture is with ϕ as the fugacity coefficient and K as the K-factor. Mikyška and Firoozabadi 16 derived a new thermodynamic function for the evaluation of the equilibrium ratio via minimization of the Helmholtz free energy that uses the specific volume, temperature, and mole fractions as its primary variables and eliminates the need for knowing the equilibrium pressure and for solving the state equation for the stable volume. They proved that the following relationship exists between the K-factor and the volume function coefficient for the liquid and vapor phases: in which Φ i is the volume function coefficient of the component i and can be computed analytically as a function of temperature, specific volume, and mole fractions via We define as the molar specific value of the volume function, such that, instead of using Eq. (7), the natural logarithm of K-factors can be calculated by It can be shown that this (molar) specific volume function is related to the fugacity coefficient via By substituting the general cubic EoS (2) for the evaluation of the partial pressure term in the integral (9), the following expression is obtained for this new thermodynamic function: where g i is The equality of chemical potentials and component material balances can be systematically expressed by means of K-factors in such a way that the vapor mole fraction θ is determined by the classic Rachford-Rice equation whereẑ i is the overall mole fractions of component i in the feed. Then, molar compositions of the liquid and vapor phases are obtained: This formulation leads to a remarkable reduction in the number of variables in isothermal flash calculations; we know the overall composition of the feed hence, by knowing the n K-factors, we can

Reduction method
The basic idea of all reduction methods is to calculate the K-factors in a lower-dimensional hyperspace spanned by parameters that are independent of the number of components in the mixture.
According to the classical theory of reduction 20 , such reduced parameters can be obtained by decomposing the symmetric matrix β ij = 1 − κ ij that represents the binary interactions into matrices composed of its eigenvectors and eigenvalues, that is, Definingŝ ki ≡ s ki √â i as entries of the reduction matrix with size of m × n, we can express g i (i = 1 . . . n) in Eq. (13) as as a function of the reduced parameters Similarly, the energy parameter a of the mixture in Eq. (3) can be calculated from these reduced parameters via Then, an equation for evaluation of the molar specific value of the volume functions, ψ i , can be derived by substituting g i (i = 1 . . . n) and a into Eq. (12) using Eqs. (18) and (20): where coefficients h are functions of q k (k = 1 . . . m), b, and v: Because the entries of reduction matrixŝ ki andb i are equal in the liquid and vapor phases, all K-factors can be computed from: We note that these h-based reduced parameters are Lagrange multipliers of the classical reduced parameters, similar to the reduced parameters introduced by Nichita and Graciaa 15 . Hence, the reduced-space iteration has a better condition number and will converge faster than other methods 22 .

Thermodynamic relations for non-isothermal flashes
For non-isothermal flash calculations, it is necessary to compute additional thermodynamic quantities such as the specific molar enthalpy, internal energy, and heat capacities at constant volume and pressure. They are typically calculated as a summation of the ideal part, which is here evaluated as a function of temperature using the 9-coefficient NASA polynomials 23 , and the excess part obtained from the state equation using the reduced parameters. Overall mixture quantities are computed through where η ∈ {u, h, c v , c p } are specific internal energy, enthalpy, and heat capacities at constant volume and pressure. The molar specific internal energy of the liquid or vapor (superscripts L and V are not repeated for brevity) is computed via where u ig i is the ideal gas (NASA polynomial) molar specific internal energy of pure component i; a is obtained from Eq. (20) and its first temperature derivative is where sgn(ϑ i ) is the sign function of variable ϑ i = 1 + c(ω i )(1 − T /T ci ) and its value is equal to plus one for ϑ i > 0 and equal to minus one otherwise. The molar specific enthalpy of the mixture is defined as where p is either known or computed via Eq. (2). The molar specific heat capacity at constant volume for a multi-component mixture can be computed via Here, c ig v,i is the ideal gas molar specific heat capacity at constant volume for the component i, which is computed as a function of temperature using NASA polynomials, and the second derivative with The molar specific heat capacity at constant pressure of the mixture is computed from the thermodynamic relation

Algorithm 1: VT and PT flash calculations Result: K-factors of a multi-component vapor-liquid equilibrium
Step 0: Estimate initial values of K-factors using the input values or via the Wilson's correlation in case of blind flashes; while convergence criteria not met do Step 1: Calculate θ by solving the Rachford-Rice equation (Eq. 14); Step 2: Determine molar compositions x and y (Eq. 15) and then compute parameters q k and b for both phases (Eqs. 19 and 4); Step 3: Compute molar specific volumes v L and v V using pressure equality and volume constraint equations in case of VT-flash and two state equations of liquid and vapor in case of PT-flash; Step 4: Evaluate Jacobian matrix and update the reduced principal variables for the Newton-Raphson iteration or update the principal variables via their definitions in case of the successive substitution method; Step 5: Update K-factors (Eq. 23) and check the convergence criteria. end where the derivatives of pressure with respect to the specific volume and temperature are and 4 Numerical algorithms

Isothermal flashes
In this section, numerical solution procedures for two important isothermal phase splitting cases, PT and VT-flashes, are presented. In Algorithm 1, we need to estimate the initial K-factors at Step 0: if there is no promising data available (blind flash), Wilson's correlation is commonly employed for the initialization of the iteration. This is straightforward if the pressure and temperature are known as in the PT-flash; in the case of a blind VT-flash, however, the pressure is unknown. In this case, one could estimate the pressure from the state equation of the mixture by using the total specific volumev, temperature T , and overall mole fractionsẑ i as an input, but this will result in negative pressures in many cases. A simple remedy is to set a minimum value in pressure estimation 19 , or to employ the initialization method based on the vapor pressures of the components 24,25 . We propose to use the geometric average of the pressures of the dew and bubble points estimated as where p sat i is the vapor pressure of the pure component i, which can be estimated from Raoult's law and Wilson's correlation. In Step 1, we need to solve the Rachford-Rice Eq. (14) to determine the vapor mole fraction.
Usually, a Newton method is coupled with a bisection method for reasons explained by Michelsen and Mollerup 21 . In order to preserve the fully quadratic convergence rate of the Newton method, where For any starting value σ 0 in the range of (0, +∞), monotonic convergence of the Newton iteration is guaranteed for one of these two functions. The estimated value of σ is updated via where G and H are derivatives of G and H with respect to σ: and The Newton iteration is repeated with σ new until the convergence criteria is met. The vapor mole fraction is then obtained via θ = (c 1 + σc n )/(1 + σ). In Step 2, molar compositions of the liquid x i and vapor y i for (i = 1 . . . n) are computed using Eq. (15). Then m reduced parameters q k and mixture co-volume parameter b are obtained using Eqs. (19) and (4 )for both phases. In Step 3, the energy parameter a is computed for the liquid and vapor phases using their m reduced parameters via Eq. (20). For the case of PT-flash, in which the value of the equilibrium pressure p is given, the specific volume is then computed for the vapor and liquid phases separately.
For general cubic EoS Eq. (2) this means to find the roots of the cubic equation that is written below for the liquid phase: where The same equation is holds for the vapor phase. We use Cardino's algorithm to determine all roots of Eq. (42). If more than one real root is found, the root associated with the lowest Gibbs free energy is selected 21 .
For the VT-flash, in which the value of the total molar specific volumev is given, we first compute the molar specific volume of one phase from the volume constraint (1 − θ)v L + θv V =v and then substitute it into the pressure equality equation. The resulting equation is a quintic function of the other phase specific volume that is given below for the liquid phase: where Here, parameters α i (i = 1 . . . 5) are computed via the following expressions using the liquid and vapor co-volume and energy parameters: Since there is no analytical solution, Eq. (44) has to be solved by iterative methods to obtain v L . We use a Newton method with a starting point very close to the co-volume of the mixture in this study. Afterwards, the vapor's specific volume is obtained through the volume constraint and the associated Jacobian matrix are calculated, in which δ αβ is the Kronecker delta function. Next, the resulting set of linear equations J∆ h ∆ = e, can be solved by using the Gauss elimination method with partial pivoting to In order to find the analytical expressions of the entries of the Jacobian matrix (48), we used the classical m + 2 reduced parameters including q k (k = 1 . . . m), b, and θ as the helping variables in the derivative chain rule for the required partial derivatives The required partial derivatives of the coefficients h are obtained via Eq. (22). The derivatives with respect to the reduced variable q k are In addition, the derivatives with respect to the co-volume of the phase are where a is computed via Eq. (20) as a function of reduced parameters The derivatives respect to the specific volume of the phase are Next, the partial derivatives of the specific volume in Eq. (50) are obtained through the implicit function theorem. For PT-flashes, we directly utilize the general cubic EoS (2) for each phase as follows: ∂v j /∂q j k = −(∂p/∂q k ) j /(∂p/∂v) j , for j = L, V with ∂p/∂v from Eq. (33). By using the relationship between the a and q k , we can compute ∂p/∂q k as Moreover, the derivatives with respect to the co-volume of the mixture are with It is obvious that, in PT-flashes, the partial derivatives of the specific volumes with respect to the vapor mole fraction are zero, i.e.
For utilizing the implicit function theorem for the VT-flashes, we define the function f ≡ p L − p V and then compute the required derivatives as where partial derivatives of f can be computed using the chain rule. For instance, when j = L, we obtain ∂f /∂v L = (∂p/∂v) along with ∂v V /∂v L = ∂q V k /∂q L k = ∂b V /∂b L = (θ − 1)/θ. Subsequently, the partial derivatives of specific volumes respect to the vapor mole fraction are computed through where ∂f /∂θ for the liquid and vapor phases are Finally, partial derivatives of the reduced parameters q k (k = 1 . . . m) as well as b with respect to principal variables h ∆ β (β = 1 . . . m + 2) can be obtained via their definitions: for all m reduced parameters and for the co-volume parameter In both equations, we need the derivatives of the phase mole fractions z i , which is equal to x i and y i for liquid and vapor phases, with respect to the principal variables. Using the Rachford-Rice equation and the definition of the equilibrium ratio, we obtain and The partial derivative with respect to principal variables is expressed as follows for all K-values: and for the vapor mole fraction and the index the index in the range from β = 1 to m + 2: In Step 5, the logarithm of equilibrium ratios are computed from the updated principal variables via Eq. 23 and the following convergence criterion is checked: We propose and use ε k = 10 −2 for the initial SSI and ε k = 10 −10 for NRI, but one SSI step is usually enough for the most of cases. If the solution is not converged, we jump back to step 1 with the new K-values.

Non-isothermal flashes
In this section, the numerical solution method for the HP-and UV-flashes are explained. The main idea is to use the most appropriate isothermal flash (that is, PT for HP and VT for UV) and iterate its input temperature in such a way that the specific internal energy (UV) or enthalpy (HP) converge to the given value. The numerical procedure is summarized in Algorithm 2 and explained in more detailed in the following: In Step 0, the temperature of the mixture is estimated. To provide an initial guess at regions close Algorithm 2: UV and HP flash calculations Result: Equilibrium temperature Step 0: Estimate the initial value of temperature; while convergence criteria not met do Step 1: Execute one VT-flash or PT-flash according to the availability of the specific volume or pressure and the latest available temperature; Step 2: In case of UV-flash, compute the specific internal energy and c v of the mixture, or in case of HP-flash, compute the specific enthalpy and c p of the mixture; Step 3: Update the temperature and check the convergence criteria; end to the critical point or near phase boundaries, one can estimate the temperature by considering the mixture as single-phase and then iterate the EoS for the given specific internal energy or enthalpy.
In Step 1, we perform an isothermal flash calculation using the method that most closely corresponds to the targeted non-isothermal problem, that is, we perform a PT-flash in the case of the and for the case of HP-flashes: Using the line search L ensures global convergence of the algorithm and renders the temperature initial guess less important. Subsequently, the relative error is computed, i.e. ε r = |(ĥ − h mix )/ĥ| or ε r = |(û − u mix )/û| for HP or UV-flashes, respectively. Steps 1 to 3 are repeated until the   Table 3: Molar composition of vapor and liquid in equilibrium at the states selected for the detailed analysis of the convergence of the flash algorithms. convergence criterion is satisfied, for the calculations presented in this paper until ε r < 10 −10 .

Numerical results
We have developed a Fortran implementation of the proposed flash algorithms for the four discussed isothermal and non-isothermal flash calculations, and tested it for a large number of different multicomponent mixtures and different cubic EoS. The selected representative cases that we will present discuss in the following use the PR EoS and the values for the critical temperatures, critical pressures and acentric factors that are listed in Table 1.

Convergence behavior and robustness
Two mixtures with specified compositions including a synthetic condensate gas and synthetic oil are  Next, VT-flashes have been conducted along isochores drawn on the phase diagram. Selected isochoric lines are drawn in Fig. 1   The overall thermodynamic properties at these points are listed in Table 2 and results for the molar composition of the vapor and liquid in equilibrium are shown in Table 3. These values are equal for all types of flash calculations.  Figure 3: PT-flash convergence for the Y8 and MY10 mixtures at the points marked in Fig. 1.    Figure 6: UV-flash convergence for the Y8 and MY10 mixtures at the points marked in Fig. 1.
iterations, with only two iterations difference between points close to and far from the extreme conditions. Figure 5 shows  First, PT-flash calculations were carried out for a mixture with 26.54% ethane and we record the total computational time 100 × 100 states in the pressure-temperature range that is enclosed by the black box in Fig. 7. Then, the mixture internal energy and specific volume that were computed by the PT-flashes are used for executing the corresponding UV-flashes. In order to assess the performance of the proposed UV-flash at conditions that similar to what we typically encounter in CFD simulations, initial guesses of pressure and temperature were computed by adding random perturbations to the true values, that is:

Computational time
where r is a random number generated in the range [−0.5, 0.5]. The perturbation amplitudes ∆T and ∆p are set to 20 [K] and 20 kPa, which corresponds to the maximum change that we can expect between two subsequent time steps in CFD simulations.
The results are shown in Fig. 8. The computational time for the current PT-flash algorithm is always lower than the highly optimized reference method. The difference becomes more significant as the number of components is increased, which shows the importance of reduction methods for the both iso-thermal and non-isothermal flashes. Surprisingly, we also measure a performance gain for the two-component mixture, where the number of variables is not reduced by the new method.
In this case, the reduction method acts as a preconditioner and reduces the number of required iterations for the PT-flash. Furthermore, it should be noted that the computational performance the UV-flash based on the VT-flash is much less sensitive to the amplitude of the imposed pressure perturbation ∆p than the conventional method based on the PT-flash. For instance, the conventional method becomes more than five times slower for ∆p = 400 kPa, whereas the overall time needed for the new method remains unchanged.