A Gibbs free energy minimization based model for liquid – liquid equilibrium calculation of a system containing oil, brine, and surfactant

. Accurate and reliable phase equilibrium calculations of microemulsion systems are of great importance. This study deals with the thermodynamic modeling of Liquid – Liquid Equilibrium (LLE) of a system including oil ( n -decane), brine (containing CaCl 2 salt), and ionic surfactant (sodium dodecyl sulfonate). Two models of UNIQUAC and UNIQUAC + Debye – Hückel were used for thermodynamic calculations. The LLE experimental data were utilized to estimate the binary interaction parameters of UNIQUAC model and the adjustable parameter, b , of the Debye – Hückel model. The thermodynamic model calculates the microemulsion phase ’ s compositions by minimizing the Gibbs free energy of the LLE system using a combination of genetic algorithm and fmincon function in order to prevent local minima. The thermodynamic modeling results show an appropriate agreement with the experimental data. Accordingly, the presented model of this study can be used as a suitable method to investigate the liquid – liquid equilibrium of systems containing oil, water, and surfactant.


Introduction
Petroleum products (oil, gas, and fuels), colloidal systems, and those which result from biomass conversion are among the fluids that have received much attention in various industries [1,2]. Because of the high demand for oil around the world, more attention is being paid to Enhanced Oil Recovery (EOR) methods [3]. Given the pressure drop in conventional reservoirs, many methods are used to recover the residual oil in reservoirs [4]. For instance, water flooding has been extensively considered as an effective method in the secondary oil recovery [5]; however, given the high capillary forces and heterogeneity of reservoir, nearly 70% original oil in place remains in the reservoir [6]. Therefore, the importance of tertiary oil recovery methods is specified. In these methods, different chemicals including surfactants, polymers, and so on along with seawater are injected to oil reservoirs [7][8][9]. To this end, different surfactants, including cationic, anionic or even a mixture of these two, are widely used [10][11][12][13]. As a result of the surfactants injection, the surfactant is initially located due to its structure between the oil and water phases, which leads to the formation of a three-phase system including aqueous brine, oil as well as, surfactant, and also reduces the interfacial tension [14][15][16][17]. Increasing the surfactants' concentration results in the formation of a continuous phase, or the microemulsion. Finally, a two-phase system consisting of oil and microemulsion phases, which are in equilibrium, are formed [17,18].
According to the previous researches, one can conclude that surfactant plays a very important role in the LLE of oil, brine, and surfactant systems. Therefore, accurate and reliable phase equilibrium calculations of microemulsion systems to facilitate the EOR conditions are of great importance so that the exact number of phases and the values of the components are determined in equilibrium [19][20][21][22]. As indicated in the other researches, there are two methods for calculating liquid-liquid equilibrium: in the first method, the isoactivity equations are solved under mass balance constraints [23,24], and in the second method, the Gibbs free energy minimization is addressed [25]. In the second method, the results are dependent on initial estimates, and poor estimates may lead to wrong solutions [26]. For instance, in 1999, Yushan and Zhihong proposed an algorithm for liquid-liquid phase equilibrium calculation. They utilized tangent plane distance function criterion for stability analysis so that this approach was considered as sufficient and necessary condition. Moreover, activity coefficient models, namely UNIQUAC and NRTL were used to formulate Gibbs free energy functions [27]. In 2018, Dadmohammadi et al. employed the NRTL, mNRTL1, and mNRTL2 activity coefficient models for the characterization of LLE systems [28]. Their study showed the efficiency of the modified NRTL models in characterization of LLE systems. Also in 2016, Li et al. unlike the use of minimizing the Gibbs free energy (coupled with stability test), solved isoactivity equations for their LLE calculations using NRTL model, and conducted their study on two-phase LLE systems. They concluded that this approach could be applied to multiphase systems [26]. It should be noted that both Equation of State (EOS) and activity coefficient models have been implemented in different researches for characterization of the phase behavior of LLE [29][30][31][32][33][34][35] and vapor-liquid equilibrium [36] systems. Riazi and Moshfeghian presented a thermodynamic approach using NRTL activity coefficient model. They used Debye-Hückel theory to calculate the long-range interactions. They investigated the LLE of a system containing water, oil, and surfactant. For the adjustable parameter of the NRTL model, the value of 0.1 was chosen. In addition, only two interaction parameters were optimized and others considered equal to zero [37].
In this study, the LLE of oil, brine, and surfactant systems was studied through Gibbs free energy minimization method on the data provided in Riazi and Moshfeghian's research [37]. To do so, two models of UNIQUAC and UNIQUAC + Debye-Hückel were used, and consequently the results were compared with the previous research.

Experimental data
As mentioned earlier, the experimental (phase equilibrium) data of a system including oil (n-decane), brine (containing CaCl 2 salt), and an ionic surfactant, named sodium dodecyl sulfonate, were obtained from Riazi and Moshfeghian's research [37]. The studied system includes two phases of oil and microemulsion, so that in each of the two phases, all three components of oil, brine, and surfactant are available. It should be noted that the salt presents in water is CaCl 2 with a concentration of 500 mg/L (Ca ++ ), which is equal to 0.138 wt%. In addition, the experiments were carried out at a temperature of 20°C and atmospheric pressure. Three concentrations of surfactants with values of 0.5 wt%, 1 wt%, and 2 wt% were used in the experiments, and in each set of experiments ten different oil volumes ranging from 5% to 90% of the total volume of system were utilized. The initial mixture compositions for three sets of experiments as well as the physical properties and mole fractions of the components in microemulsion phase are presented in Tables 1-3, respectively. It is noteworthy that the / parameter represents the molar ratio of the oil phase to the total number of mole in the system, the amount of which alters from 0 to 1. It should be noted that numbers 1-3 are considered for oil, water, and surfactant molecules, respectively.
The structural parameters for the functional groups are reported in Table 4 [38,39]. The molar excess Gibbs and activity coefficient functions of UNIQUAC model are as follows [40]: where r, q, and q 0 indicate molecular structure parameters of a pure component, expressing the volume and surface areas of molecule, respectively, where with the exception of water and some small alcohols, q = q 0 . Also, the subscript k and superscript i in equations (7) and (8) refer to the group k and the occurrence of the group k in molecule i, respectively. The parameter Z, which implies the V 20 , d 20 , and e are the molar volume at 20°C, liquid density at 20°C, and the dielectric constant, respectively. coordination number, is assumed as a fixed value of 10. The variable È* represents the segment fraction. In addition, h and h 0 denote the area fractions of a simple and a hydrogen bonded molecule, respectively.

Debye-Hückel theory
The ionic activity coefficient for electrolyte is as follows [41]: where Z + and Z À are the positive and negative charges of the electrolyte. The parameter b is an adjustable parameter, and I is the ionic strength which is defined as follows: where m i and Z i stand for the molality of ion i (mol/kg of water) and the valances of ions, respectively. In Debye-Hückel theory, the A c coefficient is calculated as follows: In the calculations, it is assumed that the surfactant is fully dissociated into positive and negative ions just in the microemulsion phase, and the dissociation does not occur in the oil phase. The experimental data of this study were measured at 20°C [37]. The constant values utilized in equation (12) are presented in Table 5. The calculated value for A c is equal to 1.1667.

Gibbs free energy minimization
Calculations were performed based on the Gibbs free energy minimization [19]. Total Gibbs free energy of a mixture of n c components and n p phases is obtained from the following equations: where, where n k , n k i , F, and z i are the total number of moles in phase k, the number of moles of component i phase k, the total number of moles in input feed, and the molar fraction of component i in the input feed, respectively. In the case of using the Gibbs energy functions for the LLE calculations, the following equation is utilized to calculate the Gibbs energy of the mixture: where the superscript, L(k) denotes to liquid phase k, and c i represents the activity coefficient. Through the minimization of the above equation, equilibrium concentrations are calculated in each phase. In the above equation, the mole number of each component in each phase is unknown that can be calculated by minimizing. To achieve this goal, change of variables, b k i (for i = 1, 2, . . ., n c ; k = 1, 2, . . ., n p À1), is used which is as follows: The unknown variables are b k i , which are calculated by minimizing. The following limitation must be applied to avoid negative or imaginary answers:

Liquid-liquid equilibrium calculation
In this study, to find the accurate values of b k i in the minimization of Gibbs energy and the binary interaction parameters of UNIQUAC and UNIQUAC + Debye-Hückel models, a combination of genetic algorithm and fmincon function (fmincon function within MATLAB environment) was utilized which leads to calculation of logical answers in all points.
One of the methods to solve constrained and unconstrained optimization problems is to use genetic algorithm based on natural selection, a process that leads to biological evolution. The way it works is that the genetic algorithm frequently modifies a population of individual solutions. At each stage, the genetic algorithm randomly chooses individuals from the current population as parents and also it utilizes them to produce children for the next generation. Consequently, throughout consecutive generations, the population "evolves" toward an optimal solution. Genetic algorithms can be applied to solve different optimization problems that cannot be implemented using standard optimization algorithms. These include problems where the objective function is stochastic, discontinuous, nonlinear, or non-differentiable [42,43]. In addition, the fmincon function is used to find a constrained minimum of a scalar function whose variables start from an initial estimation.
As we know, in the calculations of the liquid-liquid equilibrium, the equilibrium condition is as follows: where O and ME denote to the oil and microemulsion phases, respectively. In equation (21), the activity coefficients are calculated based on the mole fraction. The problem with this is that the Debye-Hückel model is based on molality. Therefore, the conversion of the unit of activity coefficient from molality to mole fraction is required. To do so, the equation provided by Panah et al. can be used [44], which is as follows: where x, c (x) , m, c (m) are the mole fraction, mole fraction based activity coefficient, molality, and molality based activity coefficient, respectively. In the above equation, the term of "1000/18.01528" is the mole number of water contained in 1 kg. When using the UNIQUAC + Debye-Hückel model for the surfactant, the UNIQUAC model is used to calculate the activity in the oil phase and the Debye-Hückel model is applied for the microemulsion phase. Therefore, the normalization of the Debye-Hückel model is necessary and it is proposed as follows: ln where c Ã 3 and c 1 3 are the calculated activity coefficient of the surfactant by the Debye-Hückel model and the surfactant activity coefficient at infinite dilution, respectively. In the calculations, the value of c 1 3 was obtained by fitting which is reported in Table 8.
To find the binary interaction parameters of UNIQUAC activity coefficient model, the adjustable parameter of Debye-Hückel model (b), as well as the surfactant activity coefficient at infinite dilution (c 1 3 ), equation (21) should be solved. Afterwards the Gibbs free energy minimization will be performed.
In UNIQUAC + Debye-Hückel model, the activity coefficient of n-decane and water in both phases are calculated using UNIQUAC equations. The activity coefficient of the surfactant in the oil phase is also calculated using UNIQUAC equations, but the activity coefficient in the microemulsion phase is calculated using Debye-Hückel equations.
In the last step of the present study, Absolute Average Deviation (AAD) is used to find the quality of the predicted data compared to the experimental data. The AAD% is defined as below: where N is the number of tie lines and X is the variable whose error is calculated. In addition, the superscripts "exp" and "cal" refer to experimental and calculated data, respectively.

Results and discussion
The structural parameters for oil, water, and surfactant molecules in UNIQUAC activity coefficient model are calculated using equations (7) and (8). These parameters are presented in Table 6. As Teh and Rangaiah mentioned in their study, in order to have precise phase equilibrium calculations, an accurate global minimum is required [19]. As we know, although the genetic algorithm optimizes the parameters globally, this optimization is not accurate, since different initial guesses for parameters give different answers. Moreover, functions such as fmincon optimize the parameters locally. To solve these problems, in the present study, the combination of these two algorithms was used so that both binary interaction and b k i parameters are first optimized by genetic algorithm and then by fmincon function, this procedure leads to accurate results.
The binary interaction parameters for UNIQUAC and UNIQUAC + Debye-Hückel models are presented in Tables 7 and 8, respectively. The UNIQUAC and UNIQUAC + Debye-Hückel models have 6 and 7 binary interaction parameters, respectively. In UNIQUAC + Debye-Hückel model, full dissociation of the surfactant into two positive and negative ions in the microemulsion phase is assumed, but in UNIQUAC model, no dissociation is considered.
The calculated values of the compositions (%) in the microemulsion phase as well as / values (%) for UNIQUAC and UNIQUAC + Debye-Hückel models are presented in Table 9. In addition, for a better comparison, for example the predicted values of the present study for surfactant composition as well as / values are shown in Figures 1  and 2, respectively, as compared to the experimental data and the model presented by Riazi and Moshfeghian [37].
It should be mentioned that the mixture of water and CaCl 2 is not considered as a pseudo component. The concentration of salt in the initial mixture is very low and its presence is ignored in the oil phase.
The AAD% of the present modeling results are presented in Table 10 in comparison with the results acquired by Riazi and Moshfeghian [37].
Due to the low concentrations of surfactants, the addition of electrolyte term or, in other words, the consideration of long-range interactions, has little effect on  [37]. Moreover, the comparison between the results obtained by Riazi and Moshfeghian [37] with the results of this study shows a much better performance of the UNIQUAC and UNIQUAC + Debye-Hückel models, especially in predicting the oil composition in the microemulsion phase. It is necessary to mention that all three models have poor performances in predicting the molar fraction of the oil component. However, UNIQUAC and UNIQUAC + Debye-Hückel models predict the molar fraction of the oil component with a higher accuracy than the NRTL model used by Riazi and Moshfeghian [37]. Nevertheless, in the analysis of the obtained results, this point should be considered that the NRTL model requires four fitting parameters to describe the proposed LLE system; while the obtained lower errors by the UNIQUAC and UNIQUAC + Debye-Hückel models were originated from more fitting values  [37] offers a more successful performance in predicting the molar fraction of the surfactant component, and there is an appropriate agreement with the experimental values. The UNIQUAC and UNIQUAC + Debye-Hückel models represent an excellent performance in predicting È. However, in the NRTL model used by Riazi and Moshfeghian [37], a very high error is observed. One of the goals pursued in modeling studies could be to reduce the fitting parameters to the lowest possible number as the high number of fitting parameters can reduce the validity of the model. In 2015, Dadfar and Biria studied the LLE of the system with the same experimental data used in the present study [21]. They used 35 fitting parameters. However, the error obtained by them were not considerably better.
It should be mentioned that, as the oil composition is not reported in the literature [37], the oil is considered as n-decane. Moreover, in modeling study carried out on the same system, it is common to consider oil as n-decane [21].

Conclusion
Gibbs free energy minimization based modeling of the LLE conditions for a system containing oil, brine, and surfactant was investigated in this work. This modeling approach employs UNIQUAC and UNIQUAC + Debye-Hückel activity coefficients models. A combination of genetic algorithm and fmincon function was used which leads to calculation of reasonably accurate results. According to Table 10, one can note that the investigated model in this study displays an acceptable overall performance in predicting the molar fractions of components as compared to a previous study reported in literature so that the experimental and computational values were slightly different. It is concluded that the addition of electrolyte term has little effect on the improvement of results given the low concentrations of surfactants. In fact, considering the electrolyte term, the AAD% value for the mole fractions of the components and / changes less than 1%. The model can be used as a robust method to investigate the liquid-liquid equilibrium of systems containing oil, water, and surfactant.