Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443
Dissipative particle dynamics simulation of the phase behavior and microstructure of CTAB/octane/1-butanol/water microemulsion Zhengxia Chen, Xiuli Cheng, Haishi Cui, Ping Cheng, Hongyan Wang ∗ College of Chemistry, Jilin University, Changchun 130012, People’s Republic of China Received 10 October 2006; received in revised form 10 January 2007; accepted 14 January 2007 Available online 21 January 2007
Abstract In this paper dissipative particle dynamics (DPD) has been used to predict the phase behavior and microstructure of Cetyltrimethylammonium bromide (CTAB)/octane/1-butanol/water system. The dynamics process for the formation of water/oil (w/o) microstructure can be well reproduced by the method. The effect of octane on the microstructure has been studied by varying octane concentration. By the analysis of the simulated interface tension and microstructure, their phase diagram has been fabricated, which is consistent with experimental result. The mesoscopic properties of microemulsion, such as micelle shape, density distribution of water in the micelle, interfacial tension, can inform us the relation between micelle bulks and component ratio. By adjusting the bulk and figure of the micelles, the expected nanometer structure material can be synthesized with controlled structure according to the phase diagram. © 2007 Elsevier B.V. All rights reserved. Keywords: Dissipative particle dynamics; Microemulsion; Phase behavior; Microstructure
1. Introduction Microemulsion is transparent, isotropic and thermodynamically stable dispersion system of oil and water. The size of liquid drop in the dispersed phase is in the range of 8–80 nm. The quaternary system of surfactant, cosurfactant (a short chain linear alcohol), oil and water has many important features, so it is one of the extensively studied microemulsion systems. In the last decade, a number of experimental techniques have been used to investigate the properties and microstructures of microemulsion, such as infrared spectroscopy (IR), nuclear magnetic resonance (NMR), quasi-elastic light scattering (QELS), conductivity measure, etc. [1–4]. These techniques can give us in-depth understanding for some features of microemulsion. However, its precise knowledge of microstructure and phase behavior is still missing. The reason for this is that the complicated multicomponent mixture would bring formidable obstacles to experimental determination study of microstructure. Thus, it is necessary for us to seek one new method to study microemulsion and understand the
∗
Corresponding author. Tel.: +86 431 5168470; fax: +86 431 5175863. E-mail address: wang
[email protected] (H. Wang).
0927-7757/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.colsurfa.2007.01.022
mechanism of the phenomena taking place inside the complex system. In recent years, some computer simulation methods have been devoted to the study on surfactant/water/oil system, especially for vesicle or suspension. The vapor/liquid interface of three-dimensional Lennard–Jones fluids has been studied by molecular dynamics (MD) and Monte Carlo (MC) [5–8] methods, but it is difficult to investigate the liquid/liquid interface of Lennard–Jones fluids containing several thousands atoms, which would spend much computational time. Mesoscopic simulation techniques, such as the dynamics self-consistent-field theory (SCFT) [9–12] and dissipative particle dynamics (DPD) [13–20] have been applied to study the morphology evolution in many fluid systems. These complex fluids include surfactant monolayer, monomer mixtures, polymer solutions, diblock copolymer, etc. In these methods, the surfactant and polymer chain are commonly treated as a coarse-grained bead-spring model. Each bead represents a group of atoms. Mesoscopic simulation can treat a wider range of length and time scales compared with MD simulation. Yuan et al. [17] put forward the simple model in which sodium bis(2-ethylhexyl) sulfosuccinate (AOT) is represented by one-head and two-tail beads tied with a harmonic spring. They have drawn ternary phase diagram of (AOT)/water/iso-octane by the changes of interface
438
Z. Chen et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443
tension between water and oil phase. Qian et al. [18] have successfully applied DPD to study the phase formation of cyclic diblock copolymer c-AmBn (m + n = 20). The simulated result is consistent with the experimental result. Sarah et al. [19] have reported the phase behavior of amphiphilic polymers by DPD method. They can predict the formation process of poly(ethylene butylene)–poly(ethylene oxide) emulsion in water and methylcyclohexane, which cannot be obtained by experiment. Ryjkina et al. [20] have applied DPD simulation method to investigate the phase behavior of non-ionic surfactants. Thus, DPD is one effective simulation method to investigate the phase behavior of the complex fluids. However, the studies on the behavior and microstructure of the quaternary microemulsion by DPD are scarce. In this paper one quaternary system composed of Cetyltrimethylammonium bromide (CTAB)/octane/1-butanol/water has been researched. Firstly, we study the dynamics process of forming reverse micelles. Secondly, the effect of octane concentration on microstructure is investigated in detail. Under the fixed octane concentration, the biggest single-phase microemulsion has been checked to validate the applicability of the selected force field and models. Under the effective force field, the simulated phase diagram for the system is compared with experimental result [21]. These results can provide some interface information in the microemulsion, such as interfacial tension, bulk and figure of dispersed drop, which is very helpful for the study of nanometer particle with controlled structure in future. And it can inform us the reaction mechanism taking place in the little reactors among the complex system, which cannot be shown by experimental techniques. 2. Computational method Dissipative particle dynamics is a relatively new simulation method to study the hydrodynamics behavior of complex fluids over long length and time scales [22–24]. The complex fluids include dispersed system, emulsion and macromolecule solution, etc. The method is based on the dynamics of a set of soft beads, of given mass and size, which interact with other beads via soft potentials. Here, we describe an outline of the model and the evolution algorithm of the DPD. A bead represents a cluster of molecules and its motion is assumed to be governed by Newton’s laws. dri = vi dt mi
where aij is the repulsion parameter between particle i and j and rij = ri − rj , rij = |rij | and rˆij = rij /rij . The dissipative force FD is proportional to the relative velocities of the two beads and acts as to reduce their relative momentum. The random force FR provides an energy input into the system and builds a thermostat together with the dissipative force. They are given by the following equations: FijD = −γωD (rij )(ˆrij · vij )ˆrij
(5)
FijR = σωR (rij )θij t −1/2 rˆ
(6)
where vij = vi − vj , σ is the noise amplitude, γ the dissipative parameter, ωD (rij ) and ωR (rij ) the weighting functions which become zero for r ≥ rc and θ ij is the random Gaussian variable with zero mean. It is convenient with simulation units m = kB T = rc = 1. Herein, the surfactant (CTAB (2)) is divided into two DPD beads as shown in Fig. 1. The C16 H33 fragment and N+ (CH3 )3 fragment are treated as hydrophobic tail (T) and hydrophilic head (H) bead, respectively, which are tied together by a harmonic spring. Water (1), 1-butanol (3) and octane (4) are represented by one bead W, COS, O, respectively. There are several methods suggested in the literatures to evaluate the interaction parameters for DPD simulation. The quantitative structure-property relationship (QSPR) is based only on the solubility parameter [25]. It is not appropriate for the present system because of a partial charge in CTAB molecule. Molecular mechanics method is usually more accurate but takes a considerable amount of computational work. In compromising the accuracy and computation time we adopt the Monte Carlo method to evaluate the free energy of two species from their pair contact energies. In the first step, following Ryjkina et al. [20], the mixing energy of two corresponding fragments i and j Eijmix can be obtained by
(1)
dvi = fi dt
(2)
where ri , mi , vi and fi are position vector, mass of the particle, velocity and total force on the bead, respectively. The force between each pair of particles is made up of a conservative term FijC , a dissipative term FijD and a random term FijR . The effective force Fi acting on bead i is given by Eq. (3) Fi = FijC + FijD + FijR (3) i=j
The conservative force FijC is a soft-repulsive interaction, which is linear up to a cutoff distance rc : ⎧ ⎨ a 1 − rij rˆ , r < r ij ij ij c rc (4) FijC = ⎩ 0, rij ≥ rc
i=j
i=j
Fig. 1. DPD particle structures model of water (1), CTAB (2), 1-butanol (3) and octane (4); Cetyltrimethylammonium bromide (CTAB). “H” denotes the hydrophilic head group (N+ (CH3 )3 ) and “T” denotes the hydrophobic tail (C16 H33 ).
Z. Chen et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443 Table 1 Interaction repulsion parameters aij in simulation system (unit: kB T)
W H T COS O
W
H
T
COS
O
25.000 25.339 151.519 61.997 103.237
25.339 25.000 177.823 90.471 143.607
151.519 177.823 25.000 26.737 25.943
61.997 90.471 26.737 25.000 28.075
103.237 143.607 25.943 28.075 25.000
439
sensitive to the dissipative parameter γ in the range of 2–32kB T. Herein, the dissipative parameter γ is set to a value of 4.5kB T. All DPD beads belonging to the same molecule are connected by a loosely bounded spring with a spring force constant K = 4. This spring constant controls the stiffness of the molecule, but it is not very sensitive to the simulation result. The 20,000 DPD step is adopted in this research. 3. Results and discussion
using equation, Eijmix
3.1. Dynamics process of formation reverse micelles
1 = [Zij Eij (T ) + Zji Eji (T ) − Zii Eii (T ) 2 − Zjj Ejj (T )]
(7)
where the Zii , Zij , Zi and Zjj are the calculated values of the coordination numbers for each pair of fragments. The mean pairinteraction energy Eij (T) can be obtained from Monte Carlo calculations as shown in equation [26,27], dEij P(Eij )Eij exp(−Eij /kT ) Eij (T ) = (8) dEij P(Eij ) exp(−Eij /kB T ) P(Eij ) represents the probability distribution of pair-interaction energies. These calculated averaged mixing energies are used to obtain the Flory–Huggins parameter xij (T) as shown in equation, xij (T ) =
Eijmix
(9)
RT
The xij (T) parameter is used to determine the maximum repulsion force between the particle pair i, j at given particle density, ρ = N/V, given in equation aij (T ) = aii + 3.27xij (T )
for ρ = 3
(10)
The aii term is derived from the compressibility of pure component i (aii = 75 kB T/ρ) [22]. To calculate the mean pairinteraction energy Eij (T), the assignment of partial charge on each individual atom is required. The calculated interaction energy is very sensitive to these partial charges. Herein, the optimized geometries and partial charges of atoms of CTAB molecule are obtained by the quantum mechanical method INDO/2. The molecule is partitioned into two DPD fragments as described. The averaged interaction energy is calculated based on these fragments. The COMPASS force field is chosen to carry out Monte Carlo calculations. Later the DPD simulation confirms the suitable choices of the force field and charge assignment. COMPASS is the most recent and accurate force field for the calculation of molecular interactions in aqueous systems. This force field is specifically parameterized for the modeling of fluids [28,29]. The calculated aij parameters are given in Table 1. All computational works are performed on the SGI workstation using the software Cerius2 . The size of simulated box is set to 20 × 20 × 20 in rc units, which contains a total of 24,000 beads. The bead density of all systems is set to 3 (in relative units, related to rc ). In the literature [22], the simulation test has shown that the simulated results are not particularly
The system is composed of 7.5% water, 6.75% CTAB, 76.5% 1-butanol and 9.25% octane molecules (mol%). In Fig. 2, the change of microstructures with increasing DPD steps is shown by coarse grain model, in which the 1-butanol and octane are not shown for clarity from Fig. 2a–e. Fig. 2f displays the existent state of the every component at 20,000 steps. Obviously, at 100 steps (Fig. 2a), CTAB and water beads are distributed uniformly in the simulated cell. With the increase of DPD steps, surfactant molecules begin to absorb around the water clusters. When the DPD step is 2600 (Fig. 2b), CTAB molecules begin to form the initial reverse micelles in oil system. Then, the figures of reverse micelles are changed from out-of-order state into regular sphere when the DPD steps are increased from 8000 to 10,000 (Fig. 2c and d). Comparing Fig. 2e with Fig. 2d, the figures and the distances between two spherical micelles change little, which indicates one stable system at 20,000 DPD steps. In Fig. 2e, all blue hydrophobic tails of CTAB are pointing outside and the green hydrophilic heads are pointing inside of the reverse micelles, which indicate less repulsion force of water with head than that of water with tail. The phenomena are consistent with the repulsive parameter value in Table 1. In Fig. 2f, azury octane beads and pink 1-butanol beads as medium are well distributed in the microemulsion system. Thus, the dynamics process of forming W/O microstructure can be easily reproduced. And the microstructure at 20,000 DPD steps is selected to study the phase behavior of different system in the following research. 3.2. Effect of octane concentration on the microstructure of microemulsion The octane concentration has a significant influence on the stability of water/CTAB/octane/1-butanol system. The effect of the octane concentration on the microstructure has been investigated as follows. Fig. 3k, i, l–o are the microstructures at 1.4, 2.8, 5, 10, 13 and 20% of octane concentration. The simulation system is set 20 × 10 × 10, and the interface tension at x-direction [24] can be obtained at different octane concentration in Fig. 4. When the concentration of octane is 1.4%, the micelles are the rodlike (Fig. 3k). With increasing concentration of octane, the additive spherical micelles are formed (Fig. 3i). Nevertheless, after 5% of octane concentration, phase separation has happened between octane and water (Fig. 3l and m). While at 13% of octane concentration, sandwich liquid crystal phase is formed in Fig. 3n. With sequentially increasing octane con-
440
Z. Chen et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443
Fig. 2. Dynamics process of the forming W/O system (1-butanol and octane are not displayed). From (a–e), the aggregates are after 100, 2600, 8000, 10,000 and 20,000 DPD steps, respectively. Red, water; blue, hydrophobic beads; green, hydrophilic beads of surfactant; but the last one (f) including all particles is depicted at simulated steps 20,000; azury, octane; pink, 1-butanol. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
centration, phase separation happened again in Fig. 3o. These changes are consistent with that of interface tension in Fig. 4. With the increase of octane concentration in the range of 0–5%, the interface tension between water and octane decrease gradually, and then increase gradually resulting in phase separation. At 13% of octane concentration, one peak of surface tension is shown because of the forming liquid crystals. And the liquid crystals are changed into two phase again at 20%. Thus, the microstructure of microemulsion has been affected by the concentration of octane. And octane concentration of 2.8% (namely,
10% percent in the mixture of octane, 1-butanol and CTAB) is selected for the investigation of the biggest single-phase region in following simulation. 3.3. Microstructures of the biggest single-phase region The experimental phase diagram of CTAB/octane/1butanol/water system at 298 K is shown in Fig. 7b [21]. It shows that the microstructures of the largest single phase in quaternary system are water-in oil (W/O) (a–d), bicontinuous (B.C) (e, f,
Fig. 3. Microstructures of for the different octane concentrations: (k) 1.4%, (i) 2.8%, (l) 5%, (m) 10%, (n) 13% and (o) 20%. Blue, hydrophobic beads and green, hydrophilic beads of surfactant; azury, octane; pink, 1-butanol. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Z. Chen et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443
441
Table 2 Simulation microstructures for CTAB/1-Butanol/Octane/water system Concentration (mol %)
Microstructures
Item (dots)
Water
1-Butanol
Octane
CTAB
Simulation
Experiment
a d f g i
7.5 26 43 54 72
76.5 43.2 27.9 22.5 13.5
9.25 7.9 5.7 4.6 2.8
6.75 23.4 23.4 18.9 11.7
W/Oa B.Cb B.C B.C O/Wc
W/O W/O B.C B.C O/W
a b c
Water-in oil. Bicontinuous phase. Oil-in water.
Fig. 4. Interface tension in octane/water/CTAB/1-butanol system. The interfacial tension is relative intension of DPD units, not the true one.
g), oil-in water (O/W) (h, i) when the octane concentration is 10% in the mixture of octane, 1-butanol and CTAB. The phase behavior for the dots from a to i (Fig. 7b) has been investigated through experiment. We also studied these dots by DPD method in this section. In Fig. 5, the phase behavior and microstructures is simulated with isosurfaces model, which connects the points with equal density field values of water. The microstructures are shown in Table 2. The change of interface tension is shown in Fig. 6. From Fig. 5a, W/O micelles are formed in the oil phase when the water concentration is 7.5%. The interface tension is very low as shown in Fig. 6. From 26 to 54% (Fig. 5d, f and g), the microstructure is B.C, in which oil and water phases interpenetrate with each other in the role of surfactant. The interface tension increases quickly and then keeps one constant. Thus, one shoulder dot can be found. With the continuous increase
Fig. 6. Interface tension in octane/water/CTAB/1-butanol system with the increase of water concentration.
Fig. 5. Microstructures of the biggest single-phase region by isosurfaces model (connect points with equal density field values of water) for the different water concentrations: (a) 7.5%, (d) 26%, (f) 43%, (g) 54% and (i) 72%.
442
Z. Chen et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443
could be explained that the role of CTAB is confined in the DPD method when CTAB is separated into the head and tail as two beads. In experiment, the long hydrocarbon chain, tail, can provide the steric hindrance for water. In DPD method, the tail is only considered as one bead as big as head. In a word, it is easier and more economical to reproduce the microstructure of microemulsion by DPD method than by experiment. 4. Conclusions We have successfully applied dissipative particles dynamics to simulate the phase behavior and microstructures of microemulsion. The simulated result is in good agreement with experimental phase diagram. This paper has firstly proved that DPD method can be used in the complex quaternary system. And the information at the mesoscopic level is helpful in the application of microemulsion, which is inaccessible in experiment. Acknowledgement The authors acknowledge the National Science Foundation of China (Grant No. 50573027) for support of this research. References
Fig. 7. Phase diagram for CTAB/1-butanol/octane/water system at 298 K: (a) the simulated result and (b) the experimental result [21]. The coordinates of x, y, z denote CTAB + 10% octane, 1-butanol + 10% octane, water, respectively (percent represents the mole ratio of octane in the mixture of octane, CTAB, butanol).
of water concentration, O/W microstructure appears obviously with spherical and rodlike micelles (Fig. 5i). The interface tension increases quickly, and then decreases rapidly. Thus, we have obtained two points of phase transition, namely, one is from W/O to B.C before 26% of water, and the other is from B.C to O/W after 54%. By the analysis Table 2, the simulation microstructures are consistent with experimental results, which also validates the applicability of the selected force field and models in our paper. 3.4. Simulation phase diagram Combining the change of interface tension and microstructures, a large number of phase transition points can be confirmed by varying the ratio of CTAB/water/octane/1-butanol. By the connection of these phase transition points, the simulation phase diagram has been obtained in Fig. 7a, which is consistent with experimental result (Fig. 7b). But from Fig. 7, we can see a little difference between simulated phase diagram and experimental result in the region of B.C and W/O structure. This phenomenon
[1] L.J. Puig, J.C. S´anchez-Diaz, M. Villacampa, E. Mendiz´abal, J.E. Puig, A. Aguiar, I. Katime, Microstructured polyacrylamide hydrogels prepared via inverse microemulsion polymerization, J. Colloid Interface Sci. 235 (2001) 278. [2] M. Fanun, E. Wachtel, B. Antalek, A. Aserin, N. Garti, A study of the microstructure of four-component sucrose ester microemulsions by SAXS and NMR, Colloid Surf. A 180 (2001) 173. [3] X. Sui, Y. Chu, S. Xing, C. Liu, Synthesis of PANI/AgCl, PANI/BaSO4 and PANI/TiO2 nanocomposites in CTAB/hexanol/water reverse micelle, Mater. Lett. 58 (2004) 1255. [4] M. Shen, R. Guo, C. Chen, Phase behavior and structural properties in CTAB/c-C6H11OH/C6H5CH3/H2 system, J. Yangzhou Univ. Nat. Sci. Ed. 5 (2002) 25. [5] B. Derecskei, A. Derecskei-Kovacs, Z.A. Schelly, Atomic-level molecular modeling of AOT reverse micelles. 1. The AOT molecule in water and carbon tetrachloride, Langmuir 15 (1999) 1981. [6] L.J. Chen, Area dependence of the surface tension of a Lennard–Jones fluid from molecular dynamics simulations, J. Chem. Phys. 103 (1995) 10214. [7] M.J.P. Nijmeijer, C. Bruin, A.B. van Woerkom, A.F. Bkker, Molecular dynamics of the surface tension of a drop, J. Chem. Phys. 96 (1992) 565. [8] M.J. Haye, C. Bruin, Molecular dynamics study of the curvature correction to the surface tension, J. Chem. Phys. 100 (1994) 556. [9] J.G.E.M. Fraaije, B.A.C. van Vlimmeren, N.M. Maurits, M. Postma, O.A. Evers, C. Hoffmann, P. Altevogt, G. Goldbeck-Wood, The dynamic meanfield density functional method and its application to the mesoscopic dynamics of quenched block copolymer melts, J. Chem. Phys. 106 (1997) 4260. [10] R. Hasegawa, M. Doi, Adsorption dynamics. Extension of self-consistent field theory to dynamical problems, Macromolecules 30 (1997) 3086. [11] M. Laradji, A.C. Shi, J. Noolandi, R.C. Desai, Stability of ordered phases in diblock copolymer melts, Macromolecules 30 (1997) 3242. [12] M.W. Matsen, Cylinder↔Gyroid epitaxial transitions in complex polymeric liquids, Phys. Rev. Lett. 80 (1998) 4470. [13] R.D. Groot, Mesoscopic simulation of polymer-surfactant aggregation, Langmuir 16 (2000) 7493.
Z. Chen et al. / Colloids and Surfaces A: Physicochem. Eng. Aspects 301 (2007) 437–443 [14] C.M. Wijmans, B. Smit, Simulating tethered polymer layers in shear flow with the dissipative particle dynamics technique, Macromolecules 35 (2002) 7138. [15] D.W. Li, X.Y. Liu, Y.P. Feng, Bond-angle-potential-dependent dissipative particle dynamics simulation and lipid inverted phase, J. Phys. Chem. B 108 (2004) 11206. [16] L. Mohamed, J.A. Hore Michael, Nanospheres in phase-separating multicomponent fluids: a three-dimensional dissipative particle dynamics simulation, J. Chem. Phys. 121 (2004) 10641. [17] S.L. Yuan, Z.T. Cai, G.Y. Xu, Y.S. Jiang, Mesoscopic simulation study on phase diagram of the system oil/water/aerosol OT, Chem. Phys. Lett. 365 (2002) 347. [18] H.J. Qian, Z.Y. Lu, L.J. Chen, Z.S. Li, C.C. Sun, Computer simulation of cyclic block copolymer microphase separation, Macromolecules 38 (2005) 1395. [19] G.S. Sarah, K. Hubert, S. G¨unter, M. Christain, V. Joachim, Phase behavior of amphiphilic polymers: a dissipative particles dynamics study, Colliod Polym. Sci. 283 (2004) 284. [20] E. Ryjkina, H. Kuhn, H. Rehage, F. M¨uller, J. Peggau, Molecular dynamic computer simulations of phase behavior of non-ionic surfactants, Angew. Chem. Int. Ed. 41 (2002) 983.
443
[21] J. Xu, G.Z. Li, G.W. Zhou, S.J. Liu, Studies on phase behavior of CTAB/octane/1-butanol/water microemulsion, Det. Cosmet. 2 (2000) 35. [22] R.D. Groot, P.B. Warren, Dissipative particle dynamics: bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (1997) 4423. [23] R.D. Groot, T.J. Madden, Dynamic simulation of diblock copolymer microphase separation, J. Chem. Phys. 108 (1998) 8713. [24] C.M. Wijmans, B. Smit, R.D. Groot, Phase behavior of monomeric mixtures and polymer solutions with soft interaction potentials, J. Chem. Phys. 114 (2001) 7644. [25] J. Bicerano, Prediction of Polymer Properties, second ed., Marcel Dekker, New York, 1996. [26] M. Blanco, Molecular silverware. I. General solutions to excluded volume constrained problems, J. Comput. Chem. 12 (1991) 237. [27] C. Fan, B. Olafson, M. Blanco, S. Hau, Application of molecular simulation to derive phase diagrams of binary mixtures, Macromolecules 25 (1992) 3667. [28] H. Sun, COMPASS: an ab initio force-field optimized for condensed-phase applications-overview with details on alkane and benzene compounds, J. Phys. Chem. B 102 (1998) 7338. [29] Molecular Simulations Inc., Programme Cerius2 4.1.