Interfacial free energy of the structure II cyclopentane hydrate in cyclopentane liquid from the capillary wave fluctuation method
The interfacial stiffness of structure II cyclopentane hydrate in cyclopentane liquid has been measured by the Capillary Wave Fluctuation (CWF) method for an equilibrated system with temperature of = and pressure of , close to the upper (L-Hyd-Vap-L) Quadruple-point, on the boundary of the hydrate stable region. The hydrate-II/oil interfacial stiffness results fall broadly in the range =(35–55), with an average over all simulations giving
=. Earlier work has determined the cubic harmonic expansion of the natural crystalline oscillations in structure II hydrate. Using this work, the stiffness value translates into an orientation-averaged interface free energy of =. The latter is in excellent agreement with the experimentally based estimates of = from particle-particle adhesion force measurements . This provides further evidence for the feasibility and robustness of the simulation approach to measure interfacial free energies in molecular systems.
The interfacial free energy of gas hydrate interfaces , both toward solid substrates [22, 4, 21, 7], and fluids [5, 1, 2, 3], are of crucial concern to evaluate both the likely site of hydrate initial growth [19, 20] in process streams, and as input to meso-scale growth models [28, 24]. This parameter is also of crucial importance in understanding the detailed mechanism of kinetic hydrate inhibitors and anti-agglomerants (AAs). It has the potential to elucidate if a surrounding oil phase is needed for AA action, or whether a coating of surface active components a few (tens of) molecular layers thick is sufficient to effectively screen the hydrate-hydrate adhesion forces. Recent experiments may suggest the latter . A general simulation scheme for determining the interfacial free energy of hydrate interfaces, could become part of model-based screening for potential new compounds for use in hydrate flow assurance. It could also be useful in risk assessment programs to find potential physical locations where hydrate growth is probable. Conversely, in hydrate/mineral interfaces it can be used to help evaluate stability properties of and kinetic barriers to hydrate extraction in model reservoirs.
The physics, the simulation procedure and the analysis algorithm have been described in our previous work , and only a brief summary will be given here. Molecular dynamics simulations of elongated rectangular,slab-like hydrate/fluid topologies were conducted to detect and quantify the CWFs in the interface . The interfaces have a pseudo-1D geometry with the bi-tangential length significantly shorter than the tangential one. This is done to enhance the magnitude of the CWFs . The molecular dynamics package GROMACS  was applied, with the utilization of the OPLS forcefield [17, 18] for the cyclopentane (CPN) molecules, and the compatible TIP4P/ice  rigid, non-polarizable model for the water molecules. The equations of motion were integrated by the Leap-frog half-step scheme. A Cartesian simulation box was used with periodic boundary conditions in all 3 dimensions. The TIP4P/ice model was augmented with the Particle Mesh Ewald(PME) scheme  for long range electrostatics. The parameters used were a k-space cutoff of and 4th order B-spline interpolation. This is now the preferred approach also with rigid water models , despite the fact that The TIP4P/ice model is parameterized for plain cut-off electrostatics. A single cut-off of was used for Van der Waals interactions. The water model was kept rigid by the standard SETTLE  algorithm, whereas the CPN-bonds were kept rigid during simulations by the P-LINCS  scheme, using a 4th order matrix series expansion and a single iteration. After an annealing scheme to minimize the hydrate-slab dipole moments , and a separate equilibration of the oil phase, the two phases were brought into contact. A further equilibration with the following settings were conducted on the merged hydrate/CPN(l) system:
2.1 Equilibration settings
A total of were run with a time step of , after which the energy and pressure stabilized. Neighbor searching was performed every 5 steps. The PME algorithm was used with an initial k-space grid of 320 x 15 x 512 cells. Temperature and a semi-isotropic pressure control were done with the Berendsen scheme [8, 6] with time constants, = and = respectively. Due to the presence of the interface, the tangential and normal pressures were adjusted sequentially, first the normal and then the tangential pressure to approximate the procedure in . The compressibilities of the barostat were = and = for the directions tangential and normal to the interface respectively.
2.2 Simulation settings
A total of were simulated with a time step of . Neighbor searching was performed every 10 steps. The PME algorithm was used with a k-space grid of 240 x 12 x 384 cells. Temperature coupling was done with the V-rescale algorithm , enforcing the canonical ensemble with the same time-constant as the equilibration.
The analysis was identical to our previous work  and consisted of two stages:
Restricting attention to the interface region along the normal-to-interface coordinate.
Applying a Gibbsian dividing surface criterion to the order parameter profile in the interface region, as delimited in 1.
Stage 1 is achieved by fitting the order parameter profile to a suitable interface functional, typically a hyperbolic tangent, and using the location and width parameters of the functional to delimit a specific interval along the normal(-to-interface) coordinate. This procedure is repeated for every trajectory snapshot so it is possible to correct for moving interfaces in the course of longer simulations. This scheme helps eliminate the resolution limitation that would otherwise be incurred by a rectangular spatial mesh of the order parameter in the box. The mass density was used as order parameter throughout.
The simulation configurations are detailed in table 1. The use of the rectangular slab topology creates two interfaces in the box, which we denote the upper and lower interface according to their normal(-to-interface) coordinate. A snapshot of one simulation in its initial stage is shown in fig. 1a. Energy statistics and thermodynamics in the equilibration and simulation stages are given in table 2. The results show that the average pressure is well adjusted in the equilibration stage. Starting from a particular equilibrated state, with no barostat applied, however, leads to large deviations from the set pressure, reflecting the inherent pressure fluctuations around the average. This is particularly an issue with our extremely thin (pseudo-1D) slab topology, even in a multistage pressure relaxation. This discrepancy may indicate that our approach to approximate the scheme of pressure equilibration in  is not completely satisfactory. The CFM spectra with best fit-lines for Simulation no. 3 are given in fig. 2. The resulting fit diagnostics and estimated stiffnesses are shown in table 3. The upper and lower interfaces do not show identical physical behavior, as also reflected in the results: The two simulations at = yield surface stiffnesses of = and =. The lower interfaces tend in some cases to have liquid bridge formation of cyclopentane with its own vapor. This bridge formation (fig. 1b) occurs spontaneously within in most simulations. Only the upper interfaces remain fully immersed in cyclopentane throughout. The results from the upper and lower interface are both in agreement within with the result = from particle-particle adhesion experiments , when taking into consideration the small difference in the values of stiffness and interfacial free energy. This difference is no greater than as estimated from the fit to cubic harmonics in the hydrate-II/vacuum simulations. Due to the liquid bridge formation, the most reliable results are likely to be those from the upper interfaces in Sim. no. 3 and 4.
The scatter in the results is appreciable and many of the spectra show significant noise contamination. As in our previous work on the same simulation scheme  there are significant interface-interface dynamical correlations, even though the positive and negative cross-correlations largely average out over the whole interfacial strip, as seen in table 4. (The dispersion of the cross-correlation values are larger than their average). These correlations are indicative of interface-interface interference which plays a major role in attenuating the lower harmonics seen in the spectra fig. 2a,b. This can be avoided only by extending the system size in the normal(-to-interface) direction. The systematic differences in the free energy between the upper and lower interfaces have already been noted. The reason why cyclopentane liquid bridges form only at the lower interface in our simulations is unclear. A possible reason may be that energy is not intrinsically independent of normal coordinate in the box: In the trajectories, the lower interfaces essentially coincide with the simulation box lower boundary. This energy effect could then perhaps be due to artifacts of Ewald sum electrostatics. However, we do not have sufficient data to discount the effect as being an expression of spontaneous symmetry breaking from a saddle point in the potential, rather than a fully reproducible outcome from the initial conditions. Due to the constant volume and particle number it is observed that a formed liquid bridge will likely preclude the formation of another on the opposite hydrate/oil interface.
Further visual inspection of the trajectories indicated some surface melting of the hydrate as the capillary waves rippled across the interfaces. It was generally straightforward to judge by a combination of visual inspection and potential energy change which results were reliable and which were too contaminated to be trusted.
The melting line of the CPN-hydrate/liquid simulation system is almost identical to the near-vertical L-Hyd-L experimental 3-phase line in the P,T phase diagram, starting from the second quadruple point at =(,). We subtract =, due to the fact that the TIP4P/ice water model has ice point =. It is clear that in temperatures lower than = the capillary wave magnitudes remain small, the sufficiently large ones occurring only close to the phase transition temperature. The exception to this is the lower interface of Sim. no. 1, where a large liquid bridge is formed, and therefore the effective phase transition temperature is close to = - the melting point of CPN-hydrate in vacuum/vapor. This gives a free energy somewhere between that of hydrate/vacuum and hydrate/oil on this interface.
We can deduce from the above that, in the absence of extensive liquid bridge formation, the capillary wave method can function with a temperature fixed at or just above the hydrate/fluid melting line, provided the system remains stable (in potential energy) during data collection, but it does not function well with a temperature set significantly below the melting point of the system at the imposed pressure. Depending on the temperature increase over the phase transition point, it can be questioned, whether we in such a semi- or pseudo-stable system really measure the interfacial free energy of the solid/liquid system as opposed to a liquid/liquid or adsorbed-phase/liquid system. It is a validation of the robustness of the method that despite all these difficulties, even in the presence of some liquid bridging, the average value obtained from all simulations close to the melting line at = show such a significant overlap with the best experimental results on this system .
-  Z. M. Aman. Lowering of clathrate hydrate cohesive forces by surface active carboxylic acids. Energy & Fuels, 26:5102, 2012.
-  Zachary M Aman, Erika P Brown, E Dendy Sloan, Amadeu K Sum, and Carolyn A Koh. Interfacial mechanisms governing cyclopentane clathrate hydrate adhesion/cohesion. Physical Chemistry Chemical Physics, 13(44):19796–806, 2011.
-  Zachary M. Aman, Laura E. Dieker, Guro Aspenes, Amadeu K. Sum, E. Dendy Sloan, and Carolyn A. Koh. Influence of Model Oil with Surfactants and Amphiphilic Polymers on Cyclopentane Hydrate Adhesion Forces. Energy & Fuels, 24:5441–5445, 2010.
-  Zachary M Aman, William J Leith, Giovanny A Grasso, E Dendy Sloan, Amadeu K Sum, and Carolyn A Koh. Adhesion force between cyclopentane hydrate and mineral surfaces. Langmuir, 29(50):15551–7, 2013.
-  Zachary M. Aman, Kyle Olcott, Kristopher Pfeiffer, E. Dendy Sloan, Amadeu K. Sum, and Carolyn A. Koh. Surfactant Adsorption and Interfacial Tension Investigations on Cyclopentane Hydrate. Langmuir, 29(8):2676–2682, 2013.
-  M. E. F. Apol, A. Amadei, and H. J. C. Berendsen. On the use of the quasi-Gaussian entropy theory in noncanonical ensembles. II. Prediction of density dependence of thermodynamic properties. Journal of Chemical Physics, 109(8):3017–3027, 1998.
-  G. Aspenes. Adhesion force between cyclopentane hydrates and solid surface materials. Journal of Colloid & Interface Science, 343:529, 2010.
-  HJC Berendsen, JPM Postma, WF van Gunsteren, A. Dinola, and JR Haak. Molecular-Dynamics with coupling to an external bath. Journal of Chemical Physics, 81(8):3684–3690, 1984.
-  Giovanni Bussi, Davide Donadio, and Michele Parrinello. Canonical sampling through velocity rescaling. Journal of Chemical Physics, 126(1):014101, 2007.
-  Danxu Du, Hao Zhang, and David J. Srolovitz. Properties and determination of the interface stiffness. Acta Materialia, 55(2):467–471, 2007.
-  U. Essmann. A smooth particle mesh Ewald method. Journal of Chemical Physics, 103(19):8577–8593, 1995.
-  D. Frenkel. Simulations: The dark side. The European Physical Journal Plus, 128(1):1–21, 2013.
-  M. P. Gelfand and M. E. Fisher. Finite Size-effects in fluid interfaces. Physica A, 166(1):1–74, 1990.
-  B. Hess. P-LINCS: A parallel linear constraint solver for molecular simulation. Journal of Chemical Theory and Computation, 4(1):116–122, 2008.
-  Berk Hess, Carsten Kutzner, David van der Spoel, and Erik Lindahl. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. Journal of Chemical Theory and Computation, 4(3):435–447, 2008.
-  J. J. Hoyt, Z. T. Trautt, and M. Upmanyu. Fluctuations in molecular dynamics simulations. Mathematics and Computers in Simulation, 80(7):1382–1392, 2010.
-  W. L. Jorgensen, J. D. Madura, and C. J. Swenson. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. Journal of the American Chemical Society, 106(22):6638–6646, 1984.
-  W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. Journal of the American Chemical Society, 118(45):11225–11236, 1996.
-  Dimo Kashchiev and Abbas Firoozabadi. Driving force for crystallization of gas hydrates. Journal of Crystal Growth, 241(1–2):220–230, 2002.
-  Dimo Kashchiev and Abbas Firoozabadi. Induction time in crystallization of gas hydrates. Journal of Crystal Growth, 250(3–4):499–515, 2003.
-  Asheesh Kumar, Tushar Sakpal, Praveen Linga, and Rajnish Kumar. Influence of contact medium and surfactants on carbon dioxide clathrate hydrate kinetics. Energy & Fuels, 105(0):664–671, 2013.
-  Bjorn Kvamme, Tatiana Kuznetsova, Bjornar Jensen, Sigvat Stensholt, Jordan Bauman, Sara Sjoblom, and Kim Nes Lervik. Consequences of CO2 solubility for hydrate formation from carbon dioxide containing water and other impurities. Phys. Chem. Chem. Phys., 16:8623–8638, 2014.
-  S. Miyamoto and PA Kollmann. SETTLE—an analytical version of the shake and rattle algorithm for rigid water models. Journal of Computation Chemistry, 13(8):952–962, 1992.
-  T. Pusztai, G. Tegze, G. I. Toth, L. Kornyei, G. Bansel, Z. Y. Fan, and L. Granasy. Phase-field approach to polycrystalline solidification including heterogeneous and homogeneous nucleation. Journal of Physics - Condensed Matter, 20(40), 2008.
-  Bjørn Steen Saethre, Alex C. Hoffmann, and David van der Spoel. Interfacial free energy determination of the gas hydrate/vacuum and gas hydrate/water interfaces via capillary wave fluctuations. AIP Journal of Chemical Physics - submitted, 2016.
-  Bjørn Steen Sæthre, David van der Spoel, and Alex C. Hoffmann. Free Energy of Separation of Structure II Clathrate Hydrate in Water and a Light Oil. The Journal of Physical Chemistry B, 116(20):5933–5940, 2012.
-  Minwei Sun and Abbas Firoozabadi. New surfactant for hydrate anti-agglomeration in hydrocarbon flowlines and seabed oil capture. Journal of Colloid & Interface Science, 402:312–9, 2013.
-  Gyorgy Tegze, Tamas Pusztai, Gyula Toth, Laszlo Granasy, Atle Svandal, Trygve Buanes, Tatyana Kuznetsova, and Bjorn Kvamme. Multiscale approach to CO[sub 2] hydrate formation in aqueous solution: Phase field theory and molecular dynamics. Nucleation and growth. Journal of Chemical Physics, 124(23):234710, 2006.
-  C Vega, J L F Abascal, E Sanz, L G MacDowell, and C McBride. Can simple models describe the phase diagram of water? Journal of Physics: Condensed Matter, 17(45):S3283–S3288, 2005.
-  Carlos Vega and Jose L. F. Abascal. Simulating water with rigid non-polarizable models: a general perspective. Phys. Chem. Chem. Phys., 13:19663–19688, 2011.
-  R. L. C. Vink, J. Horbach, and K. Binder. Capillary waves in a colloid-polymer interface. Journal of Chemical Physics, 122(13):134905, 2005.
-  Luis E. Zerpa, Jean-Louis Salager, Carolyn A. Koh, E. Dendy Sloan, and Amadeu K. Sum. Surface Chemistry and Gas Hydrates in Flow Assurance. Industrial & Engineering Chemistry Research, 50(1):188–197, 2011.
|Size (after equilibration)|
|Sim. no.||HydII UC||N(CPN)||Lx[nm]||Ly[nm]||Lz[nm]|
|Sim. no.||T||Int.||# Pts.||Slope||Adj. R|||