Modified Density Gradient Theory for Surfactant Molecules Applied

Subscriber access provided by Washington University | Libraries

Thermodynamics, Transport, and Fluid Mechanics

A modified density gradient theory for surfactant molecules – applied to oil/water interfaces Xiaoqun Mu, Shun Xi, Faruk O. Alpak, and Walter G Chapman Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b00164 • Publication Date (Web): 09 Apr 2018 Downloaded from http://pubs.acs.org on April 9, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A modified density gradient theory for surfactant molecules – applied to oil/water interfaces† Xiaoqun Mu,‡ Shun Xi,‡ Faruk O. Alpak ,¶ and Walter G. Chapman∗,‡ ‡Department of Chemical and Biomolecular Engineering, Rice University, Houston, Texas ¶Shell International Exploration and Production Inc., Houston, Texas E-mail: [email protected]

Abstract

For more than a century, density gradient theory (DGT) has been developed and applied for interfacial property calculations of pure and mixed fluid systems. However, due to the local density approximation, DGT has not been applicable to amphiphilic molecules. By developing a modified DGT model with the chain contribution to the free energy, this paper extends the application of DGT model to heteronuclear chain molecules, such as surfactants. The chain contribution term is derived based on the work from iSAFT model. With the help of the Stabilized Density Gradient Theory (SDGT) algorithm developed in our previous work and the PC-SAFT equation of state (EoS), the modified DGT model is tested in water/oil/surfactant mixture systems, the results of which have been qualitatively verified with other theories and experimental data. †

Submitted to I&EC research

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction Surfactants are one of the most versatile products of the chemical industry. 1 The com-

mercial and industrial applications of surfactants range from automobiles to personal care to enhanced oil recovery (EOR). For example in upstream operations of the oil industry, crude oil droplets are trapped in the pores of rock, and much oil is left in the reservoir after water flooding. For oil droplets to deform and pass through the pores, low interfacial tension (IFT) between oil and water is required. 2 In the EOR process, surfactant along with water are injected into the oil reservoir, which is called the "surfactant flooding". The injected surfactant molecules accumulate in the interface and significantly lower the IFT between the oil and water phases. Therefore the surfactant flooding mobilizes the trapped oil in the reservoir rock pores, and effectively increases the oil recovery efficiency. The surface active nature of the surfactant molecule is due to its amphiphilic structure, which consists of a hydrophilic group and a hydrophobic group. Based on the chemical properties of the hydrophilic groups, surfactants can be classified into two major categories: the nonionic surfactant and the ionic surfactant. The ionic surfactant further breaks down into three buckets: anionic, cationic and zwitterionic surfactants. 1 In practice, different types of surfactants are often mixed to take advantage of each surfactant’s strength to achieve specific goals. 3–7 Besides the charge types on the hydrophilic group, other characteristics such as the length of the surfactant molecule and the system conditions (temperature, pressure, salinity, etc.) also have significant impacts on the interfacial behavior of surfactants. Therefore, knowledge of how surfactant behavior depends on molecular structure and system conditions is important to the design of a successful surfactant flooding process. Due to the sophistication and expense of experiments, 8,9 many theoretical and simulation methods have been developed to model the surfactant behavior in mixtures. Molecular Dynamics (MD) 10–12 and the Density Functional Theory (DFT) 13–17 are among the most

2


Page 2 of 38

Page 3 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


successful methods for surfactant modeling due to their capabilities of describing the special amphiphilic structure of surfactant molecules. MD simulation has the advantage of providing potentially atomistic details at the expense of long computational times. DFT requires a less detailed molecular model that is consistent with bulk equations of state while producing interfacial properties at a fraction of the computational time compared with MD simulations, particularly for systems with simple symmetry. One form of DFT is DGT, 18 which is famous for its simple functional form and accurate interfacial property predictions. Systematic studies using DGT to calculate interfacial properties of many systems, such as alkanes, 19–22 carbon dioxide, 23,24 nitrogen 25 and water 22,26–28 have been carried out over decades. However, there have not been many studies using the DGT model for surfactant systems. In our previous work, 22 we have tried to qualitatively model surfactant behavior in a 2-methoxyethanol/water/hexane mixture using the conventional DGT model. Although an enrichment of the surfactant molecule at the interface is observed in the equilibrium density profiles, the chemical potential of the surfactant molecule was estimated based on its concentration at each position z as shown in Fig. 1a, which is equivalent to surfactant molecule orienting parallel to the interface. This explains one limitation of the conventional DGT model: DGT simply treats every molecule in the system as occupying a single point in space regardless of its structure, and it approximates the free energy using the component’s local density and the density gradient. 29 This single-point model works fairly well for those nearly spherical molecules or homonuclear chain molecules, such as alkanes, nitrogen, water etc. as previously mentioned. The surfactant molecule with heterogeneous chain structure, however, would span multiple positions in the interface as shown in Fig. 1b. Much of the physics are missing when using the conventional DGT model for this kind of system. To overcome this limitation, a modified DGT model is developed and tested for surfactant molecules in this work. This paper is organized as follows: In Section 2, the chain 3



contribution to the free energy as an extension to the DGT model is derived based on the iSAFT model. The equation of state as well as model related parameters are introduced in Section 3. Section 4 presents modified DGT equations for a surfactant system, and the numerical methods used in the solution process. The model testing results are presented and discussed in Section 5. Detailed derivations of the chain chemical potential are included in the Appendix. 5

x 10

4

5 water hexane surfactant

3

2

1

0 0

x 10

4

water hexane surfactant

4 Density (mol/m3)

4 Density (mol/m3)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 38

3

2

1

1

2

3

z

4

Distance (m)

5

0 0

6

x 10

−9

(a)

1

2

3

z

4

Distance (m)

5

6

x 10

−9

(b)

Figure 1: In the conventional DGT model, the surfactant molecule is modeled as occupying a single point in space as shown in (a). The surfactant molecule actually spans multiple positions with a heterogeneous chain structure formed by hydrophilic segments (black circles) and hydrophobic segments (gray circles) as shown in (b). Density profiles in this example are calculated using the conventional DGT model at 323.15 K, 1.01 bar for a water/hexane/2methoxyethanol system.

2

Modified DGT model for surfactant molecules

2.1

Free energy functional for heterogeneous chain molecules

The density gradient theory was first proposed by Van der Waals, 30 and it was later reformulated by Cahn and Hilliard, 18 since when DGT has been widely used to predict interfacial properties of pure and mixed systems. In the conventional DGT model, the 4


Page 5 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Helmholtz free energy functional includes two contributions: the homogeneous Helmholtz free energy from a bulk equation of state, and the inhomogeneous part accounted by the density gradient term:

A

segment

Z " [ρ] = V

# N X 1 a0 (ρ) + vij ∇ρi · ∇ρj dV 2 i,j=1

(1)

where a0 is the homogeneous Helmholtz free energy, N indicates the number of components in the system, and vij is the influence parameter matrix. In order to extend the DGT model to describe molecules with heterogeneous chain structure, an inhomogeneous chain free energy term is required. To accomplish this, we turn to the SAFT 31–33 formalism which describes molecules as chains of spherical segments. The chain formation contribution added to the current Helmholtz free energy expression produces:

A[ρ] = Asegment [ρ] + Achain [ρ]

(2)

Eqn. 2 is a general form of the modified DGT free energy functional for heterogeneous chain molecules. One option of the chain formation energy term Achain is to use modified iSAFT of Jain et al., 34 Bymaster and Chapman. 35 In iSAFT, the chain term is a natural extension of the associating free energy functional in the limit of complete bonding. The associating free energy functional can capture the detailed structure of associating molecules at an interface. 36 According to iSAFT, the contribution to the inhomogeneous chemical potential of segment β due to the chain formation is (Eqn. 21 in Jain et al. 34 ): Z N γγ 0 X 1 XX δ ln ycontact [ρseg (r)] 1 δAchain β seg ργ (r) = ln XA (r β ) − dr kT δρβ 2 γ=1 γ 0 δρseg β (r β ) (β)

(3)

A∈Γ

In the first term of Eqn. 3, Γ (β) is the set of all bonding sites on segment β. XAβ denotes

5



the fraction of segment β that does not form a chain bonding to site A. We will take the limit XAβ → 0 for covalent bonding sites. This part is essential to account for chain connectivity between chemically distinctive segments and to enforce stoichiometry. The second term is similar to the chain contribution in orginal SAFT, 31–33 which takes into account the work necessary to bring two different chemical segments into contact. The first summation is over all the bonding segments, and the second summation is over all 0

γγ segments γ 0 that are bonded to segment γ. ycontact is the cavity correlation function for

the inhomogeneous hard spheres at contact. Its form for a bulk fluid is known from the Carnahan-Starling equation of state. 37 The inhomogeneous version of this function is not known in a tractable form and here we follow the approach in the development of the iSAFT model: 34 the cavity function is approximated by a simple geometric average of its bulk forms at different positions in the interface (refer to Eqn. 33 in the Appendix). The chain formation free energy term in iSAFT has been applied to model branched chain connectivity of star polymers and dendrimers. 35 In this work, we use it for the twogroup surfactant molecule: the amphiphilic surfactant molecule is modeled by having a head group (H) and a tail group (T ) bonded together as is illustrated in Fig. 2. Although we illustrate the surfactant molecule as having one head group and one tail group, each group is polyatomic with multiple segments: the head group is formed by water-like segments, and the tail group is formed by oil-like segments. The head and tail groups, while sitting on different positions in space, are connected by a rigid chain bond to form the surfactant molecule (H1 T1 ). The length of the rigid chain bond is σ12 = (σH + σT )/2, in which σH and σT are diameters of segments in the head and tail groups respectively. When added into a water/oil mixture, the surfactant molecules tend to partition to the interface with certain average orientation as shown in Fig. 3: the surfactant head group would preferentially be in the water rich phase, while the tail group would preferentially be in the oil rich phase. Based on the current surfactant model, Eqn. 3 is expanded and 6


Page 6 of 38

Page 7 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


water

hexane

head

tail α β

association sites

bonding sites

Figure 2: The two-group model of the surfactant molecule: the surfactant is formed by bonding a water-like head group and a hexane-like tail group together. reorganized within DGT framework. Here we define the surfactant head group as component 1 and the tail group as component 2. The local chain contribution to the chemical potential of the β group on the surfactant molecule (where α is the other group) derived from Eqn. 3 is: (r β ) µchain 1 β 12 12 = −KC − ln[ycontact (r 1 )ycontact (r 2 )] kT 2 N X µres vαj d2 ρj µα,bulk 0,α (r α ) − (r α ) − + 2 kT kT dz kT j=1 2 Z 12 1X δ ln ycontact [ρseg (r)] seg − ργ (r) dr 2 γ=1 δρseg 1 (r 1 ) 2 Z 12 1X δ ln ycontact [ρseg (r)] − ρseg (r) dr γ 2 γ=1 δρseg 2 (r 2 )

(4)

in which, KC is a constant whose definition can be found in the Appendix. Since the head and tail groups on the surfactant are bonded, the chain contribution for the head depends on the chemical potential of the tail group and vice versa. µres 0,α (r α ) is the homogeneous residual chemical potential of segment α evaluated by the equation of state with fluid densities at position r α .

7



Page 8 of 38

Non-surfactant components such as water (3) and oil (4) also have a chemical potential contribution due to the form of surfactant chains, given by: 2 µchain (r β ) 1X β =− kT 2 γ=1

Z

ρseg γ (r)

12 [ρseg (r)] δ ln ycontact dr δρseg β (r β )

(5)

Water rich phase

σ12

Interface region

θ0

Oil rich phase Figure 3: The surfactant molecules present in the interface with average orientation θ0 . Detailed derivations leading to Eqn. 4 and Eqn. 5 can be found in the Appendix. These equations give the chain contributions to the chemical potential due to the presence of surfactant molecules. Adding these equations back to Eqn. 2 leads to the modified DGT model for surfactant systems.

2.2

Chain gradient term in the bulk phase

In the bulk phase, the overall chemical potential of component β is defined as:

res chain µβ,bulk = µid β,bulk + µβ,bulk + µβ,bulk

(6)

res chain in which, µid β,bulk + µβ,bulk are given by an equation of state. µβ,bulk for surfactant head and

tail is obtained by taking Eqn. 4 to the bulk limit: µchain 1 β,bulk =− kT 2

2 2 12 µid 1 X X seg ∂ ln ybulk,contact α,bulk 12 KC + ln ybulk,contact + + ργ,bulk kT 2 i=1 γ=1 ∂ρseg i

8


! (7)

Page 9 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


And for non-surfactant component, the bulk chain chemical potential is: 2 12 µchain 1 X seg ∂ ln ybulk,contact β,bulk =− ργ,bulk kT 2 γ=1 ∂ρseg β

(8)

The bulk Helmholtz free energy due to the chain formation is: 2 µid achain 1 X seg γ,bulk 0 12 = ργ,bulk (1 − KC − ln ybulk,contact − ) kT 2 γ=1 kT

3

(9)

Equation of State and Parameters

3.1

PC-SAFT EoS

In this work, the PC-SAFT EoS is employed to provide the homogeneous free energy functions for the modified DGT model. The SAFT type of equation of state was originally developed by Chapman et al. 31–33 as an extension of Wertheim’s thermodynamic perturbation theory of first order (TPT1) to mixtures of polyatomic molecules with any number of association sites. The PC-SAFT EoS is one of the most successful SAFT versions, which was developed by Gross and Sadowski 38,39 after improving the contribution of long range attractions (dispersion) in SAFT. The main advantages of SAFT type of EoS are their capability of describing the associating interactions between molecules such as water and surfactant molecules forming hydrogen bonds, as well as describing molecular size through the chain term. The homogeneous Helmholtz free energy given by the PC-SAFT EoS is:

disp hc A0 = Aideal + Ahs + Aassoc 0 0 + A 0 + A0 0

9


(10)


Page 10 of 38

hc where Aideal is the ideal gas contribution known from thermodynamics; Ahs 0 0 and A0 are the

Helmholtz free energy due to the hard spheres and the formation of hard chains respectively; accounts Adisp is the Helmholtz free energy of long range dispersion attraction, and Aassoc 0 0 for the associating energy between molecules. For each component, there are five PC-SAFT parameters: segment number m, segment diameter σ, segment energy parameter /k, association energy Ai Bi /k, and association volume κAi Bi . For an associating component, the number of association sites Ns needs to be defined. In the test system, the PC-SAFT parameters of water and hexane are taken from the literature. 38,40 For surfactant molecule, same parameters (except for Ns ) as water are assigned to the head group, and same parameters as hexane are assigned to the tail group as shown in Table. 1. Table 1: PC-SAFT parameters of the model system. Comp.

Mi [g/mol]

mi [1]

σi [Å]

i /kB [K]

κAi Bi [1]

Ai Bi /kB [K]

Ns [1]

Head group Tail group Water n-Hexane

18.015 86.177 18.015 86.177

2.0996 3.0576 2.0996 3.0576

2.2551 3.7983 2.2551 3.7983

133.98 236.77 133.98 236.77

0.21156 0 0.21156 0

1811.4 0 1811.4 0

2 0 4 0

As mentioned, the number of association sites Ns are different for water and surfactant head group: a water molecule has 4 association sites (2 electron acceptors and 2 electron donors), and the surfactant head group only has 2 association sites (2 electron donors). In other words, water molecules can self-associate as well as associate with surfactant head groups, while the surfactant head groups can only associate with water molecules.

3.2

Binary interaction parameter kij

The bulk properties of the water/oil/surfactant mixture are needed as boundary conditions in DGT calculations. The phase behavior predictions of water/oil are challenging due to differences in interactions between water molecules and between hydrocarbon molecules. 10


Page 11 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In order to better describe the mutual solubility of water and hexane, a binary interaction parameter kij is used to mix the energy parameters in PC-SAFT EoS: √ ij = (1 − kij ) i j

(11)

For same or similar molecules, kij = 0 is a good estimate. For water and hexane mixture, however, a non-zero kij value is needed to reproduce the mixture’s mutual solubility. The optimized value kij = 0.21 for water/hexane mixture is obtained by fitting with the solubility data of hexane in water. The performance of the PC-SAFT EoS with kij = 0.21 in phase behavior prediction is tested against experimental data 41–45 in Fig. 4. With a temperature independent kij , the model prediction results match with both water and oil solubilities data pretty well, except at low temperature. Asthagiri et al. 46 and Ballal et al. 47 have also recently explained that matching mutual solubilities of water and alkanes with a constant kij is difficult without sophisticated polarization and electrostatic interactions. The kij matrix used in the calculation is presented in Table. 2. Table 2: kij matrix for the model system. kij Head group Tail group Water n-Hexane

3.3

Head group

Tail group

Water

n-Hexane

0.00 0.21 0.00 0.21

0.21 0.00 0.21 0.00

0.00 0.21 0.00 0.21

0.21 0.00 0.21 0.00

Influence parameter vi

The influence parameter vi is a unique parameter in DGT, which determines the interface width as well as the IFT of the system. Far away from the critical temperature, the temperature dependency of the influence parameter is negligible within the SAFT framework according to our experience and other literature. 48 The influence parameter values for water and hexane are obtained by fitting with the experimentally measured IFT of each pure 11



−1

10

Mutual solubility [molar fraction]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 38

−2

10

x(hexane) [exp.] x(hexane) [calc.] x(water) [exp.] x(water) [calc.]

−3

10

Water solubility in hexane

−4

Hexane solubility in water

10

−5

10

−6

10

−7

10

250

300

350

Temperature [K]

400

450

Figure 4: The mutual solubilities of water and hexane are calculated by PC-SAFT EoS with kij = 0.21. Results are compared with experimental data at different temperatures. 41–45 component at room temperature. 49,50 As for the surfactant molecule, the head group shares the same influence parameter as water, and the tail group has the same influence parameter as hexane. The vi values used in thie work are listed in Table. 3. In the mixture, the cross terms of the influence parameter vij are calculated using the following mixing rule: √ vij = (1 − βij ) vi vj

(12)

where βij is the mixing parameter for influence parameters. In our previous work, 22 βij =0 generated accurate IFT results for alkane mixtures. However for water/oil mixture, a non-zero βij value is required to reproduce the IFT data of the binary mixture. Several groups have discussed the selection of βij for water/oil mixtures: Miqueu et al. 27 reported βij =0.75 for water/methane mixture using SAFT-VR EoS, and Niño-Amézquita et al. 26 used βij =0.55 for water/methane with PC-SAFT EoS. In this work,

12


Page 13 of 38

βij =0.456 for unlike molecules is obtained after fitting with the IFT data of water/hexane mixture at 283.15 K and 1.01 bar. 51 The βij matrix for the target system is listed in Table. 4 Before adding surfactant, the DGT+PC-SAFT model with the kij and vij parameter sets are tested in a water/hexane binary mixture system. At different temperatures, the calculation results (black dash line) are in good agreement with the experimental data (red diamond) 51 as shown in Fig. 5.

52

IFT data SDGT calculation 51

IFT [mN/m]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


50

49

48

47 280

290

300

310

320

Temperature [K]

Figure 5: The IFT of water/hexane mixture is calculated by DGT (black dash line). The calculation results are compared with the experimental data (red diamond) 51 at different temperatures.

4

SDGT algorithm for modified DGT

4.1

Model equations

With the modified DGT model, we can write down the grand potential energy Ω[ρ] of the system. By minimizing the grand potential energy in the grand canonical system, the

13



Page 14 of 38

system reaches an equilibrium state: δΩ[ρ] δ = δρi δρi

A[ρ] −

Z X N

! ρi µi,bulk dV

=0

(13)

i=1

Table 3: Influence parameter vi (PC-SAFT).

vi · 1020 [J · m5 /mol2 ]

Head group

Tail group

Water

n-Hexane

1.5421

35.575

1.5421

35.575

Table 4: βij matrix for the model system. βij Head group Tail group Water n-Hexane

Head group

Tail group

Water

n-Hexane

0.000 0.456 0.000 0.456

0.456 0.000 0.456 0.000

0.000 0.456 0.000 0.456

0.456 0.000 0.456 0.000

The minimization leads to the Euler-Lagrange equation of the modified DGT in a planar interface as: µ0,i (ρ) + µchain (ρ) − i

N X

vij

j=1

d2 ρ j − µi,bulk = 0 dz 2

(14)

Solving Eqn. 14 subject to the boundary conditions: ρi (0) = ρi,A , ρi (D) = ρi,B produces the density profiles of each component in the interface, where ρi,A and ρi,B are bulk densities of component i in phase A and B respectively. D is the domain size of the calculation. This boundary value problem (BVP) with an added chain chemical potential term could not be solved using the reference fluid algorithm, which is so far the most widely used algorithm for DGT equations. 29 By introducing an evolution term, 22,52 the stabilized density gradient theory (SDGT) algorithm was developed as a general and robust algorithm to solve the DGT type of equations. Using the SDGT algorithm, Eqn. 14 becomes: N

X d2 ρ j ∂ρi + µ0,i (ρ) + µchain (ρ) − vij 2 − µi,bulk = 0 i ∂s dz j=1 14


(15)

Page 15 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for i = 1, . . . , N . ∂ρi /∂s is the evolution term such that the original Euler-Lagrange equations become time-dependent partial differential equations with s having the unit of J · m3 /mol2 . These equations are subjected to the boundary conditions:ρi (s, 0) = ρi,A , ρi (s, D) = ρi,B for all time points s. In the grand canonical ensemble, the chemical potential is conserved instead of the mass. Therefore the initial conditions of the BVP functions can be set flexibly. In this work, linear density distributions of each component across the interface are used as initial conditions:

ρi (0, z) = ρi,A +

ρi,B − ρi,A z D

i = 1, . . . , N.

(16)

Since the density profiles of the head group and the tail group are symmetric with a shift of z0 = σ12 cos θ0 in space, one differential equation for either head or tail group needs to be solved instead of two. In the following illustration, we pick the differential equation for the head group (1) as an example.

4.2

Time discretization

With proper initial conditions and boundary conditions, the BVPs (15) are solved using a time marching scheme until a stationary state is reached. To solve the SDGT equations more efficiently, a convex–concave splitting scheme is applied to the free energy function. The splitting of PC-SAFT chemical potential terms was justified in our last work 22 as:

hs hc µconvex = µid 0 0 + µ0 + µ0 hc µconvex,res = µhs 0 0 + µ0

µconcave = µdisp + µassoc 0 0 0

15



Page 16 of 38

The convex chemical potential µconvex is treated time-implicitly, while the concave chemical 0 is treated time-explicitly. The convex-concave splitting scheme not only potential µconcave 0 saved computational time, but also made the time discretization unconditionally stable with no restriction on time step size. It can, at the same time, ensure a monotonic dissipation of the free energy with respect to time. Let ∆sn denote the time step size from time level n to n + 1. Discretizing Eqn. 15 for head group in time while applying the convex–concave splitting yields: ρn+1 − ρn1 1 [ρn+1 (z)] + µconcave [ρn (z)] (z) + µconvex 0,1 0,1 ∆sn + µchain [ρn+1 (z + z0 ), ρn (z + z0 ), ρn (z)] − µ1,bulk − 1

N X j=1

v1j

d2 ρn+1 j (z) = 0 dz 2

(17)

where the chain chemical potential µchain [ρn+1 (z + z0 ), ρn (z + z0 ), ρn (z)] is given by: 1 µ1 chain [ρn+1 (z + z0 ), ρn (z + z0 ), ρn (z)] 1 = − kT KC − kT ln[y 12 (z)y 12 (z + z0 )] + µconvex,res [ρn+1 (z + z0 )] + µconcave [ρn (z + z0 )] 0,2 0,2 2 2 Z N 12 X X d2 ρn+1 [ρn (r)] δ ln ycontact 1 j n ργ (r) (z + z0 ) − kT dr − µ2,bulk − v2j dz 2 2 δρn2 (r 2 ) γ=1 j=1 2 Z 12 X [ρn (r)] δ ln ycontact 1 ρnγ (r) − kT dr 2 δρn1 (r 1 ) γ=1 (18)

In equations above, z represents the position of the head group in the interface, while z + z0 represents the position of the tail group. For non-surfactant component: ρn+1 − ρnβ β (z) + µconvex [ρn+1 (z)] + µconcave [ρn (z)] 0,β 0,β ∆sn 2 Z N 12 X X 1 δ ln ycontact [ρn (r)] d2 ρ j n − kT ργ (r) dr − µ − v (z) = 0 β,bulk βj n 2 2 δρ (r ) dz β β γ=1 j=1

16


(19)

Page 17 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.3

Space discretization

The finite difference method is used to discretize Eqn. 17 and Eqn. 19 in space. By interpolating M nodes on the space, equidistant nodes 0 =: z0 < z1 < . . . < zM +1 := D are placed on a grid with mesh size ∆z := D/(M + 1), cf. Figure 6. The number of grids that the surfactant chain bond spans is m = z0 /∆z. The second order derivative of density to distance is approximated by a central difference quotient:

d2 ρ i ρi,k−1 − 2ρi,k + ρi,k+1 (zk ) ≈ 2 dz ∆z 2

(20)

Applying the space discretization to Eqn. 17, we have for every time step n + 1: n ρn+1 1,k − ρ1,k + µconvex [ρn+1 ] + µconcave [ρnk ] 0,1 0,1 k ∆sn N n+1 n+1 X ρn+1 j,k−1 − 2ρj,k + ρj,k+1 chain n+1 n n =0 + µ1 [ρk+m , ρk+m , ρk ] − µ1,bulk − v1j (∆z)2 j=1

(21)

where n n µ1 chain [ρn+1 k+m , ρk+m , ρk ]

1 12 concave n = − kT KC − kT ln[yk12 yk+m ] − µ2,bulk + µconvex,res [ρn+1 [ρk+m ] 0,2 k+m ] + µ0,2 2 2 Z N n+1 n+1 12 X X ρn+1 1 δ ln ycontact [ρn (r)] j,k+m−1 − 2ρj,k+m + ρj,k+m+1 n − kT ρ (r) dr − v2j γ (∆z)2 2 δρn1,k γ=1 j=1 2 Z 12 X 1 δ ln ycontact [ρn (r)] − kT ρnγ (r) dr n 2 δρ 2,k+m γ=1

17


(22)


Page 18 of 38

And for non-surfactant component β, Eqn. 19 with space discretization becomes: n ρn+1 β,k − ρβ,k [ρn+1 ] + µconcave [ρnk ] (z) + µconvex 0,β 0,β k ∆sn (23) N 2 Z n+1 n+1 12 X X ρn+1 δ ln ycontact [ρn (r)] 1 j,k−1 − 2ρj,k + ρj,k+1 n − µβ,bulk − vβj dr = 0 − kT ργ (r) n 2 (∆z) 2 δρ β,k γ=1 j=1

subjected to the initial condition ρ1β,k = ρβ (0, zk ) for k = 0, . . . , M + 1.

ρi,A 0

ρi,k 1

2

3

k-1

k

ρi,B M-1

k+1

M

M+1

Figure 6: The space discretization of the domain. Due to the symmetric structure of the surfactant molecule, the head group density at position M − m + 1 < k ≤ M + 1 is equal to the density of its connected tail group, which locates in the bulk phase as shown in Fig. 7. Therefore, the densities of the head group in this area are constants:

ρ1,k = ρ2,bulk B

for k = M − m + 2, . . . , M + 1

(24)

H

0

1

2

...

k

M-m+1

T

M+1

k+m

Figure 7: The head group densities are constants in the region where M − m + 1 < k ≤ M + 1.

5

Results and Discussion In this section, the modified DGT model is implemented for interfacial property calcula-

tions of a water/hexane/surfactant mixture system. Several key characteristics of surfactant molecule are studied and discussed.

18


Page 19 of 38

5.1

Density profile calculations

One of the most important features of the surfactant molecule is the active surface partitioning, which can be reflected in the density profiles. The density profile ρi [z] of each component is obtained by solving Eqn. 21 and Eqn. 23 subject to the given boundary and initial conditions. Fig. 8 shows the solution process of the test system’s density profiles using the SDGT algorithm. 4

4

x 10

6 Head Tail Water Hexane


5 Density [mol/m3]

Density [mol/m3]

5

4

x 10

4 3 2 1

4 3 2 1

0 0

1

2

3 4 Distance [nm]

5

1

(a) s = 0.

3 2

2

3 4 Distance [nm]

5

0 0

6



3 2 1

4 3 2 1

1

2

3 4 Distance [nm]

(d) s = 220.

5

6

0 0

3 4 Distance [nm]

5

6

x 10

Head Tail Water Hexane

5 Density [mol/m3]

4

2

4

x 10

5 Density [mol/m3]

5

1

(c) s = 180.

4

x 10

0 0

4

(b) s = 150.

4

6


1

0 0

6

x 10

5 Density [mol/m3]

6

Density [mol/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4 3 2 1

1

2

3 4 Distance [nm]

5

(e) s = 250.

6

0 0

1

2

3 4 Distance [nm]

5

6

(f ) s = 280.

Figure 8: Density profiles of water/hexane/surfactant system (293.15 K, 1.01 bar) at different time steps. Linear density distributions are used as initial conditions, and the system reaches stable state after 280 time steps. In the test system, water and hexane form a liquid-liquid equilibrium system at 293.15 K and 1.01 bar, and the surfactant concentration in the aqueous phase is set to be 5e−6 (mole fraction). The setup of the system’s initial condition is shown in Fig. 8a. To make sure the boundary conditions do not change after the weighted density evaluation, bulk densities are extended by m points, which form the flat density profiles on both sides. Starting from a

19



small time step ∆s = 1e−5 (for stable convergence), the calculation is accelerated by the “time adaptation scheme”, 22 which increases the time step length based on Newton solver’s performance. It can be observed that slight changes of density profiles from Fig. 8a to Fig. 8b take a rather long time (150 time steps). And the time step length increases quickly from s = 180, which dramatically accelerates the calculation as shown from Fig. 8c to Fig. 8f. The system reaches an equilibrium state when either of the following stopping criteria P P n+1 n n is satisfied: kL2 ≤ εdiff1 , where i kρi − ρi i is over all the components, or |(Ω − P Ω n+1 )/Ω n | ≤ εdiff2 . The discrete L2 norm k · kL2 is defined as kF kL2 := ( k Fk2 )1/2 , where P k is over all the discrete grid points. The first stopping criteria means the density difference between the two time steps is smaller than εdiff1 , and the second one means the system’s grand potential energy difference between the two time steps is smaller than εdiff2 . We use εdiff1 = 1e−2 mol/m3 and εdiff2 = 1e−6 in this work. In the equilibrium density profiles (Fig. 8f), the surfactant has relatively low concentrations in both bulk phases, while it accumulates in the interface. Instead of having a single surfactant density profile as in the conventional DGT, the modified DGT generates density profiles for the surfactant head group (red line) and the tail group (green line) respectively. The peak of the head group density appears closer to the water rich phase, and the tail group’s density peak is closer to the oil rich phase. Similar surfactant density profile behavior was also reported in DFT calculations 16,17 as well as in MD simulations 53

5.2

The surfactant concentration effect

To test the surfactant concentration effect, different amounts of surfactant are added into the aqueous phase, and the equilibrium density profiles are calculated and compared in Fig. 9. From extra low (xs =1e−8) to relatively high (xs =5e−6) surfactant concentrations, the interfacial accumulation of the surfactant is enhanced as the density profile peak grows higher. This is because more surfactant molecules move to the interface as the bulk concen-

20


Page 20 of 38

Page 21 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tration increases. Based on the equilibrium density profiles, the system interfacial tension (IFT) can be calculated. IFT is the surface excess grand potential energy per unit area, which can be evaluated using the following equation: Ω − Ωbulk Asurface Z D N Z D X n(Γ (β) ) =kT ρβ (z) Dβ (z) + − 1 dz + [ares + P ] dz 2 0 β=1 0

γ=

(25)

in which, Asurface is the surface area, Dβ (z) is an integral given by Eqn. 30, n(Γ (β) ) is the total number of bonding sites on segment β, P is the system pressure, and ares = PN 1 ares 0 + i,j=1 2 vij ∇ρi ·∇ρj is the inhomogeneous residual Helmholtz free energy. The selection of the domain size D will not change the IFT calculation results, since the integration only contributes in the interface, and the surface excess grand potential energy equals to zero in the bulk region. The IFT calculation results for the mixture system at 293.15 K and 1.01 bar are depicted in Fig. 10. The system IFT decreases as more surfactant molecules are added into the system. This can be explained in combination with the density profile behaviors: more surfactant molecules move to the interface at higher surfactant bulk concentrations, and they squeeze the solvent molecules away, reducing the effective contact area for water and oil interactions while increasing water/head group and oil/tail group interactions. Therefore the IFT reduces by having more surfactant molecules in the interface. 16 At very low surfactant concentration, the modified DGT model produces the IFT that is comparable to the experimental data with no surfactant present 51 (the blue square in Fig. 10).

5.3

Tail length effect

In previous calculations, segment number m=3.0576 (same as hexane) is assigned to the tail group. To study the impact of different tail lengths, segment number m=1, 2, 3 and 4

21



4

4

x 10


4

600

Head Tail

Density [mol/m3]

500

3 2

400 300 200 100 0 0

1 0 0

1

1

2

3 4 Distance [nm]

2

5

6

3 4 Distance [nm]

x 10


5 Density [mol/m3]

5 Density [mol/m3]

6

4

4000

3000

3 2

2500 2000 1500 1000 500 0 0

1 5

0 0

6

Head Tail

3500

Density [mol/m3]

6

1

1

(a) xs = 1e−8. 6

3 4 Distance [nm]

2

5

6

3 4 Distance [nm]

5

6

4

x 10


x 10


5 Density [mol/m3]

5 4 3 2 1 0 0

2

(b) xs = 1e−7.

4

Density [mol/m3]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 38

4 3 2 1

1

2

3 4 Distance [nm]

5

0 0

6

(c) xs = 1e−6.

1

2

3 4 Distance [nm]

5

6

(d) xs = 5e−6.

Figure 9: The equilibrium density profiles of the water/hexane/surfactant system (293.15 K, 1.01 bar) at different surfactant concentrations in the water phase. are used, which are equivalent to having tail groups C1 , C3 , C6 and C9 in the surfactant. The IFT results of adding surfactants with varying tail length are displayed in Fig. 11. It can be observed that to achieve the same level of IFT reduction, a lower bulk concentration is needed for longer tail surfactant. This behavior indicates that the surfactant with longer tail length tends to be more effective in reducing IFT, and this is consistent with experimental measurements. 54,55 According to Traube’s rule, 56 the surface activity of the surfactant triples for every extra 22


Page 23 of 38

55

50

IFT [mN/m]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


45 SDGT calculation Exp. data

40

35

30 −10 10

−8

−6

10 10 Surfactant x [water phase]

−4

10

Figure 10: The IFT of the water/hexane/surfactant system as a function of the surfactant concentration in the water phase (red line) at 293.15 K and 1.01 bar. Blue square is the experimentally measured IFT of water/hexane without surfactant. 51 CH2 group. This is reflected in Fig. 11 that the four curves have roughly equal distances on the logarithmic axis. The successful description of the tail length impact makes the modified DGT model a promising tool to guide the selection of surfactant when designing different chemical processes.

6

Conclusion The application of the conventional DGT to surfactant molecules has been limited by

its approximation that no heteronuclear chain structure is considered in the model. This paper first developed a modified DGT model with an additional chain formation free energy that extends DGT to surfactant systems. The chain formation free energy, derived based on iSAFT work, serves as the free energy contribution to explain the heteronuclear nature of surfactant molecules.

23



55

50

IFT [mN/m]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 38

45

m=3

m=4

m=1

m=2

40

35

30

25 −10 10

−8

10

−6

10

−4

10

−2

10

Surfactant composition in water Figure 11: The impact of the surfactant tail length (m represents the segment number) on the surfactant effectiveness at constant temperature and pressure (293.15 K, 1.01 bar): the surfactant is more effective with a longer tail. With the help of the SDGT algorithm and the PC-SAFT EoS, the modified DGT model has been successfully applied for a water/hexane/surfactant system, where the surfactant is formed by a water-like head group and a hexane-like tail group. The density profiles as well as IFT of the surfactant system are calculated and compared with theories and experimental data at different conditions. The impact of different surfactant characteristics (such as surfactant bulk concentration, surfactant tail length) on surfactant effectiveness are studied as well. The modified DGT model enables further applications of DGT to complex systems. For the next step work, the modified DGT model will be tested in a real surfactant system with experimentally measured IFT data. Also a general and practical way of determining surfactant parameters will be studied and recommended in future studies.

24


Page 25 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgement The authors thank Shell International Exploration and Production Inc. and the Robert A. Welch Foundation (Grant C-1241) for the financial support. The authors also thank Dr. Le Wang, Dr. Florian Frank and Dr. Kai Langenbach for valuable discussions.

7

Appendix Detailed derivation of the chain formation contribution to the chemical potential in the

modified DGT model is explained in this section. We start from Eqn. 3, which is a general expression of the chemical potential due to the chain formation in DFT. Since there is only one bonding site on head (1) and tail (2) group (as is illustrated in Fig. 2), Eqn. 3 can be simplified to: 2 µchain (r β ) 1X β β = ln XA (r β ) − kT 2 γ=1

Z

ρseg γ (r)

12 [ρseg (r)] δ ln ycontact dr seg δρβ (r β )

(26)

Water (3) and oil (4) have no bonding sites, and therefore XAβ (r β ) = 1. The chain chemical potential for water and oil becomes: 2 µchain (r β ) 1X β =− kT 2 γ=1

Z

ρseg γ (r)

12 δ ln ycontact [ρseg (r)] dr δρseg β (r β )

(27)

The fraction of unbounded groups XAβ (r) of the head and tail group are calculated by: 1 I2,1 (r 1 ) exp(µ2,bulk /kT ) 1 XA2 (r 2 ) = I1,2 (r 2 ) exp(µ1,bulk /kT ) XA1 (r 1 ) =

25


(28)


Page 26 of 38

The recursive integral I2,1 and I1,2 in Eqn. 28 are given by: Z I2,1 (r 1 ) = Z I1,2 (r 2 ) =

exp (D2 (r 2 )) ∆(1,2) (r 1 , r 2 )dr 2 (29) exp (D1 (r 1 )) ∆(1,2) (r 1 , r 2 )dr 1

in which, Dβ (r β ) is: 2

1 δ(Ahs + Aatt + Aassoc ) 1 X Dβ (r β ) = − + kT δρseg 2 γ=1 β (r β )

Z

ρseg γ (r)

12 δ ln ycontact [ρseg (r)] dr seg δρβ (r β )

(30)

The first term in Eqn. 30 is the functional derivative of the inhomogeneous residual Helmholtz free energy, which is given in DGT by: X d2 ρ j δ(Ahs + Aatt + Aassoc ) res = µ (r ) − v (r β ) β βj 0,β δρseg dz 2 β (r β ) j

(31)

where µres 0,β (r β ) is the homogeneous residual chemical potential evaluated at position r β , and the second term is the gradient term. The function ∆(1,2) (r 1 , r 2 ) in Eqn. 29 defines the bonding length and orientation angle of the chain bond: ∆(1,2) (r 1 , r 2 ) = KF (1,2) (r 1 , r 2 )y (1,2) (r 1 , r 2 )

(32)

where K is a constant geometric factor which accounts for the entropic cost associated with the orientations of the two segments to form the bond. 57 y (1,2) (r 1 , r 2 ) is the inhomogeneous p cavity correlation function approximated by y (1,2) (r 1 , r 2 ) = y 12 (r 1 )y 12 (r 2 ). y 12 (r 1 ) and y 12 (r 2 ) are the homogeneous cavity correlation function for a hard sphere fluid evaluated using weighted densities. These functions are needed only at contact: 13,37

12 ycontact [ρseg (r 1 )]

1 3σ1 σ2 ς2 = + +2 1 − ς 3 σ1 + σ2 (1 − ς 3 )2

26


σ 1 σ2 σ 1 + σ2

2

(ς 2 )2 (1 − ς 3 )3

(33)

Page 27 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where ς i is given by: N π X seg ςi = ρβ (r)(σβ )i 6 β=1

(34)

13 The weighted density ρseg β (r) of segment β at position r is evaluated by:

ρseg β (r)

3 = 4π(σβ )3

Z |r−r 1 |

Modified Density Gradient Theory for Surfactant Molecules Applied

Recommend Documents