18 a 21 de novembro de 2014 , Caldas Novas-Goiás THERMODYNAMIC MODELING OF VAPOR-LIQUID EQUILIBRIUM FOR PETROLEUM FLUIDS

The phase behavior prediction is essential for the development and optimization of the hydrocarbon production in petroleum engineering applications. Well succeed projects in this field, such as management of a reservoir, multiphase flow analysis in production systems, and petroleum primary processing, depend on the correct prediction of phase behavior and the estimation of thermodynamic properties, like phase composition and compressibility factor. In this way, the present work purpose was to develop a computational tool that performs isothermal flash calculations for oil and gas multicomponent mixtures with enough precision for both project and operation. The computer program was developed in Visual Basic Application (VBA). The Peng-Robinson equation of state was used to model the vapor-liquid equilibrium. Newton-Raphson method was applied to obtain both the convergence of the Rachford-Rice equation and the cubic equation roots. The root for each phase was determined by the Gibbs energy minimization. There is a quantitative agreement between the results computed and the results provided by the commercial simulator HYSYS (AspenTech) over industrially relevant ranges of compositions, pressures and temperatures. Furthermore, the answers obtained are within precision which apply to compositional modeling for petroleum fluids. The composition and compressibility factor simulations for petroleum fluid were performed using real field conditions. The results found were investigated from the thermodynamic viewpoint and were consistent with both practice and theory.


INTRODUCTION
In general, phase equilibrium calculations involve prediction of composition, type, and number of equilibrium phases for any given system at any specified conditions of temperature, pressure and overall compositions.Phase behavior measurement and prediction are greatly important in chemical and petroleum engineering applications (Qiu et al., 2014).Phase equilibrium calculations are the foundation of a very broad range of petroleum engineering applications, such as compositional reservoir simulation, material balance models, miscibility studies, pipeline flow and separation processes (Wei et al., 2011;Gaganis and Varotsis, 2014).
Phase equilibrium calculations are perhaps the most important calculations in the petroleum industry, and equations of state (EOS) are the major thermodynamic models in these calculations (Michelsen and Mollerup, 2007).The solution of the hydrocarbon phase equilibrium problem using an equation of state requires substantial computational power due to the complexity and iterative nature of the calculations involved, particularly when the petroleum fluid needs to be described with a significant number of components (Gaganis and Varotsis, 2014).
Equations of state are basically developed for pure components, but can be applied to multicomponent systems by employing some mixing rules to determine their parameters for mixtures.The mixing rules are considered to describe the prevailing forces between molecules of different substances forming the mixture (Michelsen and Mollerup, 2007).
According to Danesh (1998), simple mixing rules, such as those that assume compounds are randomly distributed within the mixture, are quite adequate to describe hydrocarbon mixtures of petroleum fluids.
A computational tool that can predict the phase behavior using phase equilibrium calculations for hydrocarbon mixtures with sufficient precision for design applications can play a major role in the advancement of petroleum production.Cubic equations of state have been used widely to predict phase behavior, due to their performance and simplicity (Saber, 2011;Qiu et al., 2014).
In this way, this work has as purpose to develop a computational tool that performs isothermal two-phase flash calculations for multicomponent oil and gas mixtures with enough precision for both project and operation.

VAPOR-LIQUID EQUILIBRIUM MODELLING
A thermodynamic model provides the necessary relationships between thermodynamic properties and can be used in combination with fundamental relationships to generate all the properties required to perform phase equilibrium calculations (Saber, 2011).The most common thermodynamic models are the equations of state, and among them, the cubic equations of state are the most popular (Haghighi et al., 2009).

Flash Calculation Development
Flash calculation algorithms are requisite to determine the fraction of gas and liquid as well as the composition of each phase.Considering the case of two-phase equilibrium in a   -component mixture of an overall composition , where mixture mole fractions in the liquid and vapor phases are denoted  and , respectively.Also, from a thermodynamic viewpoint there is a necessary condition of equilibrium, that the chemical potential for each component be the same in both phases (Michelsen, 1982;Whitson and Brulé, 2000): Equivalently, for an isothermal system, the fugacities of the individual components have the same value, The thermodynamic model enables to calculate component fugacities, given temperature, pressure (or volume) and phase composition,  ̂ =  ̂(, , ),  ̂ =  ̂(, , ) Fugacity coefficients can be defined for the gas and liquid phases and used to define an equilibrium composition ratio.The equilibrium factors may then be used to converge upon the equilibrium split of both phases.These equilibrium factors are defined as: Note that at equilibrium, By definition of molar fraction, Introducing the liquid mole fraction () and vapor mole fraction () of the phase, the relation follows is possible, where, N is total mole number of the system.Considering the phase equilibrium calculation base,  = 1, Let the overall fraction of vapor phase be .A material balance for each component yields   relations, Substituting the Eq. ( 5) into the material balance equations yields, and, then the phase mole fractions  and  can thus be calculated from the K-factors and the phase fraction .
Finally, mole fractions in the liquid and the vapor phase must sum to unity, yielding one additional relationship, conveniently written in the form, that is the Rachford-Rice equation (Rachford and Rice, 1952),

Peng-Robinson Equation of State
The flash algorithm used in this model is based on the Peng-Robinson equation of state (PR EOS).The PR EOS was chosen here for several reasons.First, it is a simple form of the cubic equation of state, thus it is easy to implement for engineering calculations.Second, it has better performance for the prediction of vapor-liquid phase equilibrium properties over other cubic equation of state (Qiu et al., 2014).More importantly, it also has been widely and successfully used in the oil and gas industry.
The PR EOS (Peng and Robinson, 1976): Here,  is the universal gas constant,  is the molar volume, () is the temperature-dependent attractive parameter and  is the co-volume parameter.These two parameters, for the pure component , are determined applying the criteria of criticality, so that where   and   are the critic temperature and pressure of the component , respectively, and   is a temperature dependent parameter given by (Wei et al., 2011;Qiu et al., 2014), The acentric factor   for component  measures the deviation of the molecular shape from spherically symmetric structure.
The mixing rules are considered to describe the prevailing forces between molecules of different components forming the mixture.Thus, the parameters   and   of the mixture are determined by using the one-fluid type mixing rules (Whitson and Brulé, 2000), and, where   and   are overall composition of the component  and , respectively, and   is the binary interaction parameter between components  and .
By defining and Once these parameters are determined for the mixture, the compressibility factor can be obtained by solving the following equation, which is another form of the Peng-Robinson equation of state (Eq.14):

CALCULATION PROCEDURE
Successive substitution procedures have found widespread application in phase equilibrium calculation due to their robustness and simplicity.The Successive Substitution Method (SSM) is used in the flash calculator when its convergence is reasonable.The flash calculation procedure is as follows: Step 1: Feed of the input data Initially, the input data and the thermodynamic properties are supplied into the model.The required input data are: temperature, pressure and fluid composition with the molar fractions of each component.Whereas the thermodynamic properties known for each component are: molecular weight (), critic temperature and pressure (  ,   ), critic volume (  ), critic compressibility factor (  ), and acentric factor (  ).
Step 2:   -Initial estimate Initial estimates are frequently based on the additional assumption that the fluid phases form ideal solutions, in which vapor and liquid phase fugacity coefficients are dependent on pressure and temperature only (Michelsen, 1998;Haghighi et al., 2009).This implies that the K-factors are only functions  and , and the set of independent variables, reducing to Eq. 13.
The K-values in predominantly hydrocarbon systems are reasonably well approximated using Wilson's correlation (Michelsen, 1998): Step 3: Compute the vapor fraction () The Rachford-Rice equation (Eq.13) is solved using the Newton-Raphson method in the inner loop of the algorithm to find the value of the vapor fraction, .
The Newton-Raphson method is expressly (Burden and Faires, 2011): where the Rachford-Rice equation, objective function , and its derivative  ′ are defined by, The error function is computed through the following equation As criterion for convergence a tolerance 10 -5 was adopted for the error function.
Step 4: Determine the mole fractions of components in the vapor and liquid,   , and   , using the Eq. ( 10) and Eq. ( 11), respectively.
Step 5: Determine the compressibility factors by solving the cubic equation as follows:  Calculate the parameters  and  from Eq. ( 19) and Eq. ( 20), respectively, for liquid and vapor phases;  Compute the coefficients  and  applying Eq. ( 21) and Eq. ( 22), respectively, for both phases;  Substitute the values of  and  into Eq.( 23);  The compressibility factors of the liquid and vapor phases,   and   , are determined through solution of the Eq. ( 23) for each phase. Solver of the Eq. ( 23): Equation 23 can yield in either one or three real roots, but only one should be used to calculate fugacity coefficient in the next step.Therefore, Newton-Raphson method was applied for find the first root of the cubic equation.The parameters  1 ,  2 ,  3 and  4 were defined by the following coefficients: The Newton-Raphson method is expressed by: where the objective function , and its derivative  ′ are defined by The error function is computed by As criterion for convergence a tolerance of 10 -5 was considered for error function.Then, Ruffini rule was applied to obtain an equation of second order.Later, the Bhaskara formula was used to determine the other two roots of the cubic equation.Once the values of the roots are known, the root for each phase was determined from the minimization Gibbs energy () using the Eq.35 (Whitson and Brulé, 2000), The chosen root was that that yield the smaller Gibbs energy values.
Step 6: Using the calculated compressibility factors,   and   , determine the fugacity coefficient of the vapor and liquid phases using the Eq. ( 24).
Step 7: Use the SSM formula to update the equilibrium factors.Specifically: Step 8: Check the convergence condition using: The value of  is 10 -8 .The fugacity coefficients are then updated in the outer loop using the calculated mole fractions, and equilibrium conditions are checked for convergence.
The Figure 1 shows a flow chart of isothermal two-phase calculations using EOS.

FLASH CALCULATION ALGORITHM VALIDATION
In order to validate the model capability to predict the fluids phase behavior, the equations presented above and the appropriate values of binary interaction parameters between molecules have been applied to build the algorithm.The Figure 1 shows the composition of the petroleum fluid used to modeling purpose.The fluid chosen to evaluate the computational tool corresponds to light oil from a high-pressure reservoir.In the reservoir conditions, this oil presents only liquid phase.However, during the production, as the pressure decreases, two phases are formed.In this case, the knowledge of the phase behavior and thermodynamic properties calculations becomes essential.
The commercial software HYSYS (AspenTech) was applied to evaluate the model performance and accuracy and to compare the simulations performed with the aim to validate the tool developed in this work.Thus, the results obtained by Flash Calculation Algorithm (FCA) were compared with the HYSYS answers for the same fluid.The parameters and thermodynamic properties of all components are identical in HYSYS and FCA.The equilibrium constant for each petroleum fluid component was calculated using both FCA and HYSYS, as shown in the Fig. 3     In this investigation, it was also possible to evaluate the temperature influence on the fluid phase behavior.For high temperatures, it is remarkable that the vapor phase formation occurs at higher pressures.However, for all temperatures analyzed when the pressure decrease next to the surface condition (1×10 5 Pa), the vapor molar fraction values were higher than 0.8.
Figure 5 shows the compressibility factors analyze as a function of the pressure in the temperature of 358.15 K.The ideal fluid phase behavior is also present in the Fig. 5.No difference between the compressibility factors computed from the two simulators is observed.Thermodynamically, the vapor phase presents no ideal behavior for high pressures (Zv < 1), but when the pressure is next to 1×10 5 Pa, the compressibility factor is equal to one.The liquid phase shows opposite physical behavior to the vapor phase.

CONCLUSIONS
From the thermodynamic viewpoint, the results reported earlier have indicated that all the cases simulated by Flash Calculation Algorithm are consistent with the theory applied.The Peng-Robinson (PR) equation of state has been used to vapor-liquid equilibrium modeling in hydrocarbon mixtures containing CO2 and N2.The model developed predictions were validated against HYSYS simulation data over a wide range of pressure and temperature.All the above cases, a good agreement between the predictions and commercial software is observed, supporting the reliability of the developed tool.The Flash Calculation Algorithm proved to be a strong tool and very successful model for hydrocarbons multicomponent mixtures.
For future work recommendations, it is suggested that other modeling works including bubble and dew point estimation, volumetric properties and multiphase flow modeling could be built using the Flash Calculation Algorithm.

Figure 1 .
Figure 1.Flow chart of Isothermal Two-Phase Flash Calculations using equation of state.
at pressure of 1.5×10 7 Pa and temperatures of 298.15K and 348.15 K.

Figure 3 .
Figure 3. Values of Ki obtained using the FCA and HYSYS for pressure of 1.5x10 7 Pa and temperatures of 298.15K and 348.15K.

Figure 4 .
Figure 4. Comparison between the predicted results of the developed model (FCA) and simulated data using HYSYS for vapor molar fraction of the petroleum

Figure 5 .
Figure 5.The compressibility factors evaluation of the liquid and vapor phases against pressure, at T = 358.15K.
Then the fugacity of component  in phase  can be evaluated by,