**4. Computational simulation of phase behavior**

Apart from the commonly use of phase behavior laboratory test to assess microemulsion of a particular surfactant for chemical EOR application, it is viable to


**Table 1.** *Phase behavior type descriptions.*

#### **Figure 3.**

*Relationships among phase behavior laboratory test observation, microemulsion structure and Winsor's R ratio.*

use computational simulation approach to predict surfactant performance. These computational simulation approaches include the use of various molecular modeling methods such as Monte Carlo (MC) simulations, Molecular Dynamics (MD), Dissipative Particle Dynamics (DPD) and upper scale modeling methods such as Quantitative Structure-Property Relationship (QSPR) approaches. Molecular modeling tools can be used to understand microscopic effects, predict surfactants' properties and finally to optimize structures and mixtures of surfactants [19]. Molecular modeling tools, in combination with recently developed intermolecular potentials, can provide precise information about microscopic phenomena and lead to accurate estimation of thermophysical properties [20–23]. Meanwhile, QSPR is an analytical method for breaking down a molecule into a series of numerical values describing its relevant chemical and physical properties. It remains as the focus of many studies aimed at the modeling and prediction of physicochemical and biological properties of molecules [24].

#### **4.1 Quantitative structure-property relationship (QSPR)**

QSPR is an approach to relate molecular descriptors with experimental values of properties based on statistical method. Its prediction accuracy is dependent on the size and quality of database and calculation of the relevant descriptors. There are many QSPR-like terms being used for more specific situations, such as Quantitative Structure-Activity Relationships (QSAR), Quantitative Structure-Toxicity Relationships (QSTR), Quantitative Property-Property Relationships (QPPR), Quantitative Sequence-Action Model (QSAM) and Quantitative Structure-Reactivity Relationships (QSRR) [25]. QSPR models have been developed to predict properties of pure surfactants only. Development of QSPR models for mixtures of surfactants is still a challenge [26].

#### *Experimental and Computational Modeling of Microemulsion Phase Behavior DOI: http://dx.doi.org/10.5772/intechopen.101482*

Salager et al. [18] presented the use of Hydrophilic–Lipophilic Deviation (HLD) equations to attain optimum chemical EOR formulation for simple surfactant, oil and water system. The use of HLD concept to predict optimum surfactant formulation is a hybrid approach that combine HLD equations with experiments data, which is a QSPR approach. It was demonstrated that the phase behavior and optimum formulation can be manipulated with four main independent variables: brine salinity, oil alkane carbon number (ACN), surfactant parameter and temperature. However, these models are limited to simple system. Jin et al. [27] predicts the optimum surfactant salinity using HLD equation and measured parameters including the equivalent alkane carbon number (EACN), salinity, surfactant head dependent parameter, *Ksurf* value and surfactant characteristic curvature, *Cc*. The work is extended to predict the IFT behavior using the Hydrophilic Lipophilic Difference-Net Average Curvature (HLD-NAC) equation of state.

Moreau et al. [28] applied the QSPR method to predict surfactant optimal salinity based on its correlation with surfactant structures. The QSPR models have been proven in reference conditions but they cannot be extrapolated to other conditions outside the application domain.

Budhathoki et al. [29] use the HLD equation to design the ratio of surfactant mixtures to form optimal microemulsion at reservoir condition. Correct surfactant head dependent parameter, *Ksurf* and the surfactant temperature dependent parameter, *αT* values play a crucial role in the accuracy of the HLD method. Both *K* and *αT* values can be obtained via a combination use of HLD equation and a series of phase behavior laboratory work for each individual surfactant. This approach can reduce the experimental test matrix to design optimal surfactant mixture, but it is limited to surfactants with known *Ksurf* and *αT* values through extensive phase behavior laboratory work.

#### **4.2 Dissipative particle dynamics (DPD)**

Many of the interesting phenomena that occur in complex fluids occur at the mesoscale, which is roughly defined as the spatio-temporal scales ranging from 10–10<sup>4</sup> nm to 1–10<sup>6</sup> ns [30]. These scales of simulation are not feasible using MD simulation. DPD is a coarse-grained type of molecular simulation technique which could reduce the length and timescale of molecular dynamics (MD) simulations, allowing simulation of large and complex system. DPD modeling method can reach large length scales by combining molecule groups into particles or beads, and longtime scales by replacing atomistic forces with soft effective forces [31]. DPD is widely popular simulation approach due to its algorithmic simplicity and huge versatility. By varying the conservative forces between coarse-grained beads, one can readily model complex fluids such as polymers, colloids, surfactants, membranes, vesicles and phase separating fluids [30]. DPD can give insights on spatial organization of surfactants, interesting mechanistic information for films evolution or trends on surface tensions regarding structure of the adsorbed tensioactive molecules at an interface [19]. However, the challenges for a successful DPD simulation is finding robust and general methods for parameterization of the simulation system [32–33]. This is an active research area with recent approach to apply machine learning for DPD parameterization [34].

A breakthrough approach by Fraaije et al. [35] demonstrated the use of surface torque analysis in simulating surfactant phase behavior with DPD, to determine the optimal brine salinity specifically. Prior to this, QSPR statistical approach have been the only known approach for decades in determining surfactant phase behavior. Buijse et al. [36] used Fraaije et al. approach to design EOR surfactant formulation by optimizing the surfactant head and tail composition as well as the use of cosolvent. Both Fraaije et al. and Buijse et al. work are applied in simple pure oil system. Further discussion on the work of Fraaije et al. is in Section 4.2.1.

Rekvig et al. [14] varies the surfactant chain length and topology to investigate the effect of surfactant structure and composition of the monolayer on the bending rigidity. This work of Rekvig et al. is of particular interest, where the linking of bending rigidity to surfactant structure in predicting the stability of microemulsions is demonstrated. This is important because it is agreed that the bending rigidity is a key parameter in understanding structure and phase behavior of microemulsion [37]. Further details of the approach by Rekvig et al. in determining bending rigidity is discussed further in the next Section 4.2.2.

#### *4.2.1 Surface tension analysis*

The work of Fraaije et al. [35] is the only found published work for direct determination of surfactant phase behavior with theoretical foundation based on physical chemistry of microemulsion. The calculation principal is based on an observation by Helfrich [38–39] that one can calculate the surface torque density (torque) from the first moment of the molecular stress profile, provided the surface is tensionless. A positive torque implies a tendency to bend toward the oil phase and form a microemulsion with oil droplets dispersed in an aqueous phase, while a negative torque a tendency to form a microemulsion that has water droplets dispersed in an oil phase. A surface with zero torque has an indifferent tendency where the system will form the optimal or balanced microemulsion with average zero curvature. Fraaije et al. run the DPD simulation including electrostatics and ion interactions with added salt, surfactants, and oil. The relationship between mechanical coefficients and the stress profile is expressed in Eq. (2) [39–40]:

$$M\_n = \int \sigma(z) z^n dz \tag{2}$$

where *Mn* is the stress moment, *σ* is the stress tensor and *z* is the coordinate perpendicular to the surface.

In mathematics, a moment of a function is a specific quantitative measure, used in both mechanics and statistics, of the shape of a set of points. For example, for a set of data points representing mass, the *0-th* moment is the total mass, the first moment divided by the total mass is the center of mass, and the second moment is the rotational inertia. Similarly, for a set of data points representing probability density, the zeroth moment is the total probability (i.e., one), the first moment is the mean, the second moment is the variance, and the third moment is the skewness. The *n-th* moment, *Mn*, of a real-valued continuous function *f(x)* of a real variable about a value *c* is given is Eq. (3) below:

$$M\_n = \int\_{-\infty}^{\infty} (\varkappa - c)^n f(\varkappa) d\varkappa \tag{3}$$

Fraaije et al. presented a surface tension analysis namely Method of Moments to shows that torque can be calculated based on Eq. (3). This is done by calculating the torque of an microemulsion interface from the first moment, *M1*, only if the zeroth moment, *M0* (the interfacial tension) is zero exactly. Otherwise, both values of the neutral surface position and the interfacial tension has to be known. Meanwhile, in the tensionless limit, the value of *zs* is inconsequential. The simulation is run across various salinities to find the optimal salinity at which the microemulsion surface torque is zero. Note that there is no clear-cut boundary on when a positive or

*Experimental and Computational Modeling of Microemulsion Phase Behavior DOI: http://dx.doi.org/10.5772/intechopen.101482*

negative torque transitions in to a small value (around zero). Therefore, the boundaries are somewhat gradual. Fraaije et al. demonstrated that the Method of Moments is principally correct by successfully deriving known empirical coefficients of decade-old QSPR models.

It was noted in Fraaije et al. work that the torque is related to the bending rigidity and spontaneous curvature through Eq. (4):

$$
\pi \equiv kC\_0 \tag{4}
$$

where *τ* is the torque, *k* is the bending rigidity and *C0* is the spontaneous curvature. However, their approach does not allow for direct calculation of bending rigidity or the spontaneous curvature to deduce the stiffness of the interface. Furthermore, they have yet to attempt a full treatment of the compound mixture for application in actual crude oil system.

### *4.2.2 Interface fluctuation analysis*

Microemulsion structure is governed by the elastic constants, the bending rigidity, *k* and the saddle splay modulus, *ks* [41]. Bending rigidity characterizes the resistance of the interface toward bending. A low bending rigidity means large thermal undulations and low stability.

Rekvig et al. [14] used DPD to simulate surfactant monolayers on the interface between oil and water to calculate the bending rigidity by analyzing the undulation spectrum. The effect of surfactant density, chain length, adding co-surfactant and linear versus branched surfactant on bending rigidity are investigated. The results show that increase of the monolayer thickness has a larger effect on the bending rigidity than increasing the density of the layer. The bending rigidity also increases with surfactant chain length and is larger for linear than branched surfactants. Bending rigidity decrease linearly with mole fraction of short surfactants. Mixed film has a lower bending rigidity than the corresponding pure film for all mole fractions.

The work by Rekvig et al. [14] is reference with an earlier work by Goetz et al. [42] for lipid bilayer. Goetz et al. was the first to compute bending rigidity in molecular dynamics simulations. Rekvig et al. used a simple model of head, tail, water, and oil beads in DPD to capture the essential properties of ternary systems such as phase separation and adsorption. During the simulation at very low IFT, the interface is not strictly flat and undulatory waves can be observed (**Figure 4**).

These fluctuations of the interface are analyzed to compute bending rigidity. It is firstly done by characterizing the interface based on continuum theory [43]. This is then followed by the adoption of Helfrich [44] free energy of the interface [Eq. (5)]

**Figure 4.** *Undulation wave of microemulsion interface over time (T1 to T3).*

with local displacement from the average position of the interface to enable easier monitoring in simulations.

$$f\_H = \chi + \frac{1}{2}k(2H - c\_0)^2 + k\_s K \tag{5}$$

where *H* is the principal curvature of the surface at *M* (unique point), *K* is the gaussian curvature, depending on its sign, the surface will be curved like a sphere or like a saddle-splay, these number give us an idea of the surface's shape, *k* is the bending rigidity, *ks* is the splay modulus which is related to the shape of the interface (*K*) and *c0* is the spontaneous curvature.

Bending rigidity is obtained by analyzing the undulation spectrum of the interface. Based on hypotheses from Rekvig et al. [14] using equipartition principal and fast fourier transform (FFT) application to decomposes the undulation signal into different wave lengths, Eq. (5) is transformed into Eq. (6) below:

$$
\left\langle \left| \tilde{h}(q) \right|^2 \right\rangle = \frac{k\_B T}{A} \left( \frac{1}{(\chi q^2 + kq^4)} \right) \tag{6}
$$

where ~ *h* is the approximation of local displacement from the average position of the interface, q is the wave vector, <sup>2</sup>*<sup>π</sup> <sup>λ</sup>* , *λ* is the corresponding wavelength, *kB* is Boltzmann constant (1.38 � 1023 J/K), T is absolute temperature in Kelvin, *<sup>A</sup>* is the interface area and *γ* is the interfacial tension.

Given the spectral intensity, S(q):

$$\mathcal{S}(q) = \left\langle \left| \ddot{h}(q) \right|^2 \mathcal{A} \right\rangle \tag{7}$$

Combining Eqs. (6) and (7), Eq. (8) is devised:

$$\frac{1}{S(q)q^2} = \frac{\chi + kq^2}{k\_B T} \tag{8}$$

Eq. (8) is used to fit the interface undulation spectrum analysis' results to estimate the bending rigidity, *k*. In DPD simulation, all units are normalized where *kBT* is unity.

Published experiment data on bending rigidity values for microemulsion system may be used as a reference for the model predicted bending rigidity values. However, such published data is scarce or not related to monolayers. Majority of the publications are generally focused on theory. Zvelindovsky [45] mentioned that the bending rigidity for a surfactant monolayer between water and oil is usually in the range of 1–20 kBT. Martínez et al. [46] performed MD simulation for SDS surfactant in dodecane and brine system at zero salinity. The associate bending rigidity is 1.3 kBT. SDS is sodium dodecyl sulfate or sodium lauryl sulfate, sometimes written as sodium lauril sulfate. Kegel et al. [47] found that the bending rigidity for SDS surfactant with alcohol in cyclohexane and brine system is around 1 kBT. Binks et al. [48] found a bending rigidity of around 1 kBT for AOT surfactant in nonane and a brine system at optimum salinity and surfactant concentration. AOT is a twin tailed, anionic surfactant with a sulfosuccinate head group stabilized as a salt by a sodium cation. It was also reported that the bending rigidity value depends on the alkane length, where bending rigidity decreases with increase in alkane length.
