Computing Anisotropic Cavity Potential for Clathrate Hydrates - The

Subscriber access provided by ECU Libraries

A: Spectroscopy, Molecular Structure, and Quantum Chemistry

Computing Anisotropic Cavity Potential for Clathrate Hydrates Niraj Thakre, and Amiya K. Jana J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b12335 • Publication Date (Web): 14 Mar 2019 Downloaded from http://pubs.acs.org on March 15, 2019

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 28 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

The Journal of Physical Chemistry

Computing Anisotropic Cavity Potential for Clathrate Hydrates Niraj Thakre and Amiya K. Jana* Energy and Process Engineering Laboratory, Department of Chemical Engineering, Indian Institute of Technology−Kharagpur, India−721 302 Abstract Either concerning to the energy production from natural gas hydrates or employing the hydrate technology to gas handing and seawater desalination purposes, various pure and mixture clathrate hydrates need to be understood in terms of their phase behavior and stability. The ab initio methods have compelling implications in quantifying the anisotropic guest-water interactions that are responsible for the guest specific nature of the hydrates. Howbeit, the accurate cavity interactions for large guests are obscure due to computational infeasibility. With this research gap, we devise the Møller-Plesset theory with Dunning’s basis sets that are suitably advanced to the complete basis limit using the Pauling point counterpoise weight. The viability of the proposed scheme is attested with the Raman spectroscopic, second virial coefficient and viscosity data of methane and its binary hydrates with tetrahydrofuran and cyclopentane. These two promotors are chosen because they are perhaps the most useful and efficient in their class, for which the experimental cage occupancy data with methane are available to validate our new scheme of estimating the hydrate cavity potential.

Corresponding author. Tel.: +91 3222 283918; fax.: +91 3222 282250; E-mail address: [email protected] (A. K. Jana). *

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 Clathrate hydrates are naturally occurring at seafloor and permafrost region where the water and small natural gas components come in contact at moderate temperature and high pressure conditions. These hydrate reserves bear about 2×1016 cubic meters of natural gas that is perhaps twice of the other sources together.1-3 The phenomenal characteristics of the large gas uptake capacity and being specific to the guests make this compound a candidate material for the applications in gas storage4,5 and separation processes.6,7 Exploitation of enormous energy entrapped in the natural gas hydrates and their technological implications have stimulated an intensive research in fundamental understating of their phase behavior. The guest molecules stabilize the crystalline water lattice via interaction with host molecules. The formation pressures of the gases are quite high, while the hydrates of some specific low-boiling organics are stable at atmospheric conditions.8 More interestingly, the presence of a small percentage of these compounds can promote the hydrate formation by stabilizing the lattice. They enter to large cavities and allow the gases to occupy the small-sized cages and facilitate the storage at lower pressure at the expense of reduced storage capacity. This can be optimized with tuning the promoter dosage.9,10 The quantification methods for these purposes are the cage occupancy, hydration number and gas uptake parameters.11 To estimate these, the statistical thermodynamics based model of van der Waals and Platteeuw (vdWP)12 provides a robust approach that uses the Boltzmann probability factor of potential energy surface (PES) to define the spatial distribution of the guests within the cavities. Identification of the PES using Lennard-Jones (LJ) and Kihara potential using direct fit to hydrate data is debatable as it produces an infinite number of solutions.13 The second virial coefficient ( B ) or viscosity data (  ) yields a unique set of one-dimensional potential that is appropriate for point sized or spherical species. However, for asymmetric guests, the interaction energies are strongly dependent on relative angles and positions. The different energy minima are reported in the literature with relative orientations for methane14 and carbon dioxide15 interactions with water molecule. The multidimensional integral of these energies estimated from the Monte-Carlo simulations are demonstrated by Sparks et al.15,16 that outperform the LennardJones and Devonshire (LJD) approximation for determining cage occupancy. Compared to these molecular mechanics (MM) technique, the quantum-mechanical (QM) calculations offer a more


Page 2 of 28

Page 3 of 28 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


accurate estimation of energies at distance and angular configurations of the interacting molecules. It is with this intention that the present work has been undertaken. The QM methods comprise of some semi-empirical, e.g. density functional, and pure mechanistic, e.g. ab initio (second-order Møller-Plesset, MP2 and coupled-cluster CCSD(T)) methods.17,18 The remarkable thesis works on investigating the application of ab initio potentials in estimating the cavity potential are performed by Z. Cao,19-21 J. Klauda,22,23 B. J. Anderson,24,25 R. Sun26,27 and S. C. Velaga.28,29 Their works pave the way for QM techniques on a representative fragment of the hydrate followed by integrating the same for the hydrate clusters. It is noteworthy that the ab initio methods have compelling implications in estimating cavity potentials for the hydrate former gases of small size, e.g., methane, ethane, propane, carbon dioxide and argon molecules. Whereas, the same studies for higher hydrocarbons, e.g. promoters and inhibitors, are not available due to the large computation time associated with the highly accurate QM methods. With this research gap, we formulate a novel approach that estimates the cavity potential accurate up to highest available ab initio methods. It is also appealing in that the computation time required for the single point energy using our method (SCS-MP2/ACCD) is about one minute for methane hydrate and about couple of minutes for promoters (tetrahydrofuran, THF and cyclopentane, CPT) using a 64 core 2.5 GHz AMD OpteronTM processor. 2. Theory A direct estimation of the cage occupancies using these ab initio potentials offers the chemical potential difference between hydrate and empty hydrate phase,  w  H , using the vdWP theory. For equilibrium state estimation, this is equated to the chemical potential difference between pure water and empty hydrate phase,  w  L , estimated using the Holder’s equation.30 The required fluid phase nonidealities are recorded using the modified Patel-Teja (PT) equation of state31,32 and the modified universal quasichemical functional group activity coefficient33 (UNIFAC) model. The fractional occupancy ( ij ) of j th component in i th cavity is calculated using Langmuir constant ( Cij ) and the vapor phase fugacity ( f j ). The Langmuir constant can be obtained using the Boltzmann weighted cavity potential (  ) integrated over the cavity radius. The potential model of Kihara is generally used to estimate the interaction potentials of 3 ACS Paragon Plus Environment


Page 4 of 28

molecules separated at distance ( r ). The potential well depth (  ), equilibrium separation (  ) and hard core radius ( a ) are the characteristic parameters. It is noteworthy that no parameter is fitted to the hydrate data while predicting the cage occupancies and phase equilibrium conditions. 2.1 Hydrate phase model To predict the three-phase equilibrium of clathrate hydrates, the statistical thermodynamics based model proposed by van der Waals and Platteeuw12 is adopted. In this nonstoichiometric compound, the change in chemical potential by the introduction of a guest within the cavity is quantified by the partial occupation by the guests as

   w  H  RT i ln 1   ij  i j  

(1)

where, ij is the fractional occupancy of j th component in i th cavity. The parameter R stands for the universal gas constant and T is the temperature. This is calculated in eq 2 using Langmuir constant ( Cij ) described in eq 3 and the vapor phase fugacity ( f j ) of the component j . The parameter i is the number of cavities per water molecules, that has 1/23 and 3/23 values for the small and large cages, respectively. The fugacities of the component in vapor and liquid phases are calculated using the modified PT equation of state (Section 2.3). The cage occupancy can be calculated as

ij 

Cij f j

(2)

1   Cij f j j



C

4 exp   / kT  r 2 dr  kT 0

(3)

Here, the Boltzmann weighted cavity potential,  , is integrated over the cavity radius to calculate the Langmuir constant. The potential model of Kihara given in eq 4 is generally used to estimate molecular interaction potentials.


Page 5 of 28 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


 12 6   a      a    4      r  a   r a     

(4)

For the gas hydrate applications, the Kihara potential function is generally used as it considers the shape and size of the interacting species. The hydrate cavity specific formulation of Kihara potential from is as follows 

12 



6 



R r



  2 z '  11   10  a  11   5   4  a  5   R

R

r



  N  1 1  r  a

N  

R

 R 

N

R

   1  r  a  R R  

(5)



N     

(6)

Here, R and z ' are the cavity radius and coordination number, respectively. The value of coordination number corresponds to the number of water molecules per cavity that is 20 and 24 for the small and large cages of sI type gas hydrate, whereas it has the values of 20 and 28 for the small and large cages of sII type of the gas hydrate. The potential parameters  ,  and a are calculated using ab initio methods that is explained in the next Subsection 2.2. At thermodynamic state, the chemical potential difference (CPD) of water in filled and empty hydrate calculated ( w  H ) equals to the CPD of water in empty hydrate and the phase to which the hydrate is in equilibrium, i.e. liquid water ( w  L ). Holder’s equation30,31 provides a simplified estimation to estimate w  L as w  L w  L,o   RT RTo

T

P

hw dT  Vw dP  ln( x ) w w 2 RT RT To 0





(7)

The first term in eq 7, w  L,o , indicates the CPD estimated at reference point

(To  273.15( K ), Po  0(Pa)) and their values are discussed in the Section 4. The second and third term correct the quantity for temperature and pressure variations. The symbols hw and Vw denote the difference in specific heat and volume of the water in the hydrate and liquid phases. The required parameters for the calculation of these terms are presented here as, 5 ACS Paragon Plus Environment


hw  hwo 

Page 6 of 28

T

T 38.12  0.141(T  To ) dT

(8)

o

Vw  Vwo  6.695 1015 (m3mol-1Pa -1) P

(9)

The values for hwo is at standard point (To  273.15K , Po  0 Pa) are -4858.9 and -5202.2 for sI and sII type hydrates. The Vwo at standard point holds the values of 4.6 and 5.0 for sI and sII type hydrates. The last term stands for the correction for the activity of water due to the presence of foreign compounds in the liquid phase with  w as the activity coefficient and xw is the mole fraction of water in the liquid phase. The  w is calculated using the modified UNIFAC activity coefficient model that is given in Section 2.4. 2.2 Calculation of cavity potential 2.2.1 Generation of PES Grids The cell potential for gas hydrates is estimated using ab initio methods by drawing its potential energy surface for the interaction of water and guest molecule, e.g. tetrahydrofuran (THF), at its different spatial orientations within the hydrate cage. Its position in a hydrate cage can be defined by six dimensions, which include the three spherical coordinates (r ,  ,  ) of the centre of mass of THF with respect to the oxygen in the water as origin.19,28 Another three dimensions (r ',  ',  ') are the coordinates of all atoms of THF with respect to its centre of mass, among which the set of radial positions of the constituting atoms of the guest molecule with respect to its centre of mass, i.e. r ' , is kept constant to retain the geometry. In this way, the five degrees of freedoms are used to move the guest within the cavity with following geometrical constraints: the range of angles  ,   (40o , 40o ) to keep the guest molecule not too close to the cage wall and  ',  '  0  360o

to ensure the interaction with water molecules from all spherical directions.19 The r ' axis is aligned in such a way that it coincides with r axis (joining ‘O’ atom of water to the centre of mass of THF). Consequently, on rotating THF molecule along r ' axis, the dipole moments of water and THF remain in the same position. It suggests that there is no effect of this degree of freedom on cavity potential. However, with this rotation, the new orientation of H atoms of water molecule with respect to THF may have some effect on cavity potential parameters that is discussed with the results later in this article. 6 ACS Paragon Plus Environment

Page 7 of 28 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


Figure 1. Five dimensions to represent the position of THF molecule with respect to water. The ranges of the coordinates are determined by the geometrical constraints of the hydrate cavity. The grid points that are kept denser near the cage wall to capture the sharp variation in the energies at the repulsive region, construct a total of 6400 locations (Figure 1). The energies estimated at every radial distance ( Er ) are averaged to obtain the effective potential ( Eavg ) that can be fitted to the Kihara potential model. The QM methods and the basis set are chosen on the basis of the accuracy and computational feasibility of calculating the energies at 6400 orientations. With 400 points at every radial distances, the estimated energies are arithmetically averaged ( Eavg  (1 / 400) *sum( Er ) ) to obtain the effective potential that can be fitted to eq (4) for obtaining the Kihara potential parameters. 2.2.2 Ab initio Methodology for PES Generation The fundamental Hartree-Fock (HF) method to solve electronic Schrodinger equation is improved over the decades for better estimations of repulsion energies with considering spatial density of electron in the density functional theory (DFT), electronic correlations in MøllerPlesset perturbation theory (MPPT) and multi-electron wave functions in the coupled-cluster (CC) method. The MPPT and CC comprise of a good assessment of dispersion and many body interactions. The most accurate coupled-cluster method with single, double and triple level electronic excitations (CCSD(T)) has seventh power scaling of the size of the electronic system. 7 ACS Paragon Plus Environment


Page 8 of 28

Eventually, the noncovalent bond interaction energies can be accurately calculated using the MPPT of the second order (MP2) that is more accurate than DFT and close to CCSD(T) energies. As the present study deals with the larger guest molecules, the spin component scaled second order Møller-Plesset (SCS-MP2) is adopted.34 It improves the second order correlation correction ( E C ) in eq 10 by separately scaling the parallel ( E P ) and antiparallel ( E AP ) spin components as

E C  pS E P  pT E AP

(10)

where, pS and pT are the scaling parameters in having default values of 6/5 and 1/3, respectively. The QM techniques are employed using the basis sets that discretise the electronic Schrodinger equation into algebraic form. After the Pople-type basis sets, Dunning35 formulated the correlation consistent formulation as cc-pVnZ and aug-cc-pVnZ, where cc-p stands for the correlation consistent polarized, V for the valence only, nZ for the nth zeta and aug for their augmented versions with added diffuse basis functions. These show better convergence to the CBS limit as compared to the Pople-type split valence basis sets.35 The SCS-MP2/aug-cc-pVDZ method is proposed to be appropriate for the 6400 points PES generation in the present work. The energies are further corrected by analysing the discretization errors in basis sets. Amongst, the basis sets superposition error (BSSE) leads to the overestimation of the interaction energies. The counterpoise correction method is applied via calculating the energy of monomers on the dimer basis, i.e. making another molecule as ghost molecule having zero point charge. The following equations are used for estimation of raw (eq 11) and CP corrected (eq 12) interaction energies as raw  E AB  E A  E A E AB B AB A

(11)

CP  E AB  E AB  E AB E AB B AB A

(12)

where, the subscript denotes the molecule and the superscript stands for the mixed basis of the CP ) are molecules. To evaluate the BSSE, the counterpoise (CP) corrected energies ( E AB raw ) interaction energies as subtracted from the raw ( E AB


Page 9 of 28 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


CP raw BSSE  E AB  E AB

(13)

The BSSE dominates at small basis limit, whereas the effect lowers down with bigger basis sets and eventually diminishes at complete basis set (CBS) limit.20 To reach CBS, the basis sets are extrapolated using the energies recorded at the gradually increasing size of the basis sets, n. CBS These values are fitted to eq 14 for energy Eint as

Eintraw  EintCBS  A 'exp( B 'n)

(14)

The difference between the computed energy at a particular basis set and the energy at CBS is termed as the basis set convergence error CP CBS BSCE  E AB  E AB

(15)

The collective effect of these two errors (BSCE and BSSE) can be estimated using the following expression. CBS raw E AB  E AB  w  BSSE

(16)

Here, the value of Pauling point counterpoise weight ( w  1  BSCE / BSSE ) is multiplied to the BSSE and the resultant is to be added to the raw interaction energies to get the energy at complete basis. 2.3 Vapor phase model The nonideality of the vapor is addressed by the modified PT equation of state32 in terms of fugacity. This equation of state can be written in the following form,

p

RT a  v  b v (v  b )  c (v  b )

(17)

In eq 17, the parameters a , b and c have the following formulations:

a

 a (TR ) R 2Tc2 Pc

(18)

b

b RTc Pc

(19) 9 ACS Paragon Plus Environment


c

c RTc Pc

Page 10 of 28

(20)

Here, R stands for the universal gas constant. Tc , Pc and TR ( T / Tc ) are the critical temperature, critical pressure and the reduced temperature. The parameter  (TR ) is defined as    T  H2   (TR )  exp  H1 1         TC     

(21)

The values of parameters H1 and H 2 for various component can be found in the literature.33 The parameters  a , b and  c are defined as function of the parameter  c as

c  1  3 c

(22)

 a  3 c2  3(1  2 c )b  b2  1  3 c

(23)

3b  (2  3 c )b2  3 c2b   c3  0

(24)

The parameter b is the smallest positive real root of the eq 24. Rewriting the governing equation of PT-EoS, which is in cubic with respect to the compressibility factor ( z ),

z 3  (C  1) z 2  ( A  2 BC  B 2  B  C ) z  ( BC  C  A) B  0

(25)

in which, A , B and C are is computed as

A

ap R 2T 2

(26)

B

bp RT

(27)

C

cp RT

(28)

Finally, the fugacity coefficient of the pure component can be calculated as


Page 11 of 28 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


ln( )  z  1  ln( z  B) 

zM  a ln   2 RTN  z  Q 

(29)

where, the parameters M , N and Q are expressed as bc  p M  N  2  RT  (b  c) 2  N  bc  2  

(30)

1/2

(31)

bc  p Q N  2  RT

(32)

Eventually, the fugacity of the vapor phase can be calculated as

fi v  i xi p

(33)

2.4 Liquid phase model The modified UNIFAC model33 is used for activity coefficient calculation. The molar excess Gibbs free energy ( G E ) is taken as a sum of combinatorial ( GC ) and the residual part ( G R ). Eventually, for multicomponent systems, the activity coefficient can be expressed as,

 

 

ln   i   ln  iC  ln  iR

(34)

where,  i is the activity coefficient of i th component. The combinatorial part (  iC ) of the activity coefficient has the following form

 i ln  iC  ln   xi

 

 Z  i   qi ln   2  i

i    li  xi 

N

x l

j j

(35)

j 1

Here, j is the dummy index that moves over the entire component present in the system. The parameter li has the following form,

li 

Z  ri  qi   ri  1 2

(36)



Page 12 of 28

Here, Z is the lattice coordination number, which has a reasonable value of 10. The parameters i and  i are the area and segment fractions that can be evaluated as,

i 

xi qi

(37)

N

x q j

j

j 1

xi ri

i 

(38)

N

x r

j j

j 1

The pure component area ( qi ) and volume parameters ( ri ) used in the eqs 37 and 38 are expressed in the following form,

qi   vki Qk

(39)

ri   vki Rk

(40)

k

k

Here, R and Q are parameters specific to the function groups and k is the dummy index for all the functional groups present in the system. The residual part (  iR ) of the activity coefficient has the following form,

 



ln  iR   vki ln  k  ln ik k



(41)

in which,

   m km     ln  k  Qk 1  ln    m mk      m  m   n nm    n m 

(42)

Qm X m  Qn X n

(43)

n


Page 13 of 28 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


v x   v x i m i

(44)

i

Xm

i n i

i

n

The model of Kang et al.33 differs from the original UNIFAC in calculating the  mk term as, (1) (2) (3) 2  amk  amk T  amk T   T  

 mk  exp  

(1) The parameters amk ,

(45)

(2) (3) and amk for the components used in this present study are amk

available in the supplementary information of Kang et al.33 2.5 Simulation algorithm The steps for the simulation to obtain the ab initio cavity potential are performed as follows: i.

Optimize the geometry of individual molecules at the MP2/aug-cc-pVDZ level, and use these coordinates for constructing the pair of water-guest at a distance of 10 Å between their centres of mass.

ii.

Move the guest molecule with respect to oxygen at water molecule and rotate along its own centre of mass using spherical coordinates to generate a total of 6400 input files for the ab initio simulator.

iii.

Compute the single point energies (SPE) for all 6400 configurations at MP2/aug-ccpVDZ level and record the SCS-MP2 energies.

iv.

Compute the SPE of the optimized guest and water separately and record the SCS-MP2 energies to calculate the uncorrected energies using eq 11.

v.

Generate other two sets of 6400 input files written in Step (ii) by making separately the guest and water molecules as ghost molecules by simply putting a hyphen before the atomic charge of each atom and compute the SPE.

vi.

Compute the CP corrected energies using eq 12 and subsequently, compute the BSSE using eq 13.



vii.

Using the half-CP method, add half of the BSSE to the uncorrected energies to obtain the actual form of configurational energies (eq 16).

viii.

Take the arithmetic mean of configuration energies computed at same intermolecular distances (400 points) to get the 16 point energies and fit these energies at Kihara potential model eq 4 to regress the parameters. All the ab initio calculations are performed in the General Atomic and Molecular Electronic

Structure System (GAMESS-US,36 version 2018-R1-pgi-mkl) quantum chemistry software package. The thermodynamic model of hydrate is simulated in MATLAB® environment, for which a pressure search algorithm is adopted as follows, i.

Read the experimental equilibrium tempreature ( T ) for the hydrate formation and assume a corresponding guess pressure ( P ).

ii.

Using the critical properties ( Tc , Pc ) of the components present in the vapor phase, calculate the vapor phase fugacity ( fi v ) at temperature T (eq 31).

iii.

For the hydrate phase, evaluate the Langmuir constant using the Kihara potential (eq 3). Subsequently, calculate the fractional cage occupancies (eq 2) for the guest species using the Langmuir constant and the fugacity coefficient calculated in Step (ii).

iv.

Estimate the CPD for empty and filled hydrate (  w  H ) using eq 1.

v.

Calculate the CPD for empty hydrate and liquid water (  w  L ) using the Holder’s equation (eq 7) with the temperature and pressure correction terms (eqs 8 and 9) and the activity coefficient (  i ) for the last term using UNIFAC model (eq 34).

vi.

Estimate the absolute error (  w  H   w  L ) and try to minimize it by manipulating the pressure ( P ) using the secant method.

vii.

For a tolerance limit (say  w  H   w  L < 0.01), terminate the algorithm and save the value of T , P .


Page 14 of 28

Page 15 of 28 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 validate the proposed scheme of calculating the cage occupancy using a novel ab initio methodology and vapor-hydrate-liquid equilibrium (VHLE) of single and mixed hydrates, the experimental Raman spectroscopic data of Sum et al.37 and Moryama et al.38 are used for CH4 and THF+CH4 hydrates, respectively. 3. Results and Discussion The interaction energy of methane-water molecules separated at an equilibrium distance is calculated using the Dunning-type correlation consistent basis sets (Figure 2). Interestingly, we observe a pattern in estimated binding energy on increasing the size of the basis sets. Moreover, both the CP corrected and uncorrected energies are observed to progress with opposite slope towards the CBS limit. The extrapolation of the basis sets using eq 14 evaluates the CBS limit, CBS =5.36 kJ/mol, for methane-water interaction. E AB

Figure 2. Variation of interaction energies calculated using the SCS-MP2 method applied with the family of aug-cc-pVnZ basis sets. The image in inset shows the methane-water dimer configuration, for which the interaction energies are estimated. Szczȩśniak et al.14 examined the stability for different configurations of methane-water dimer and reported the relative order of interaction energies as F-H < E-H < V-O < E-O < F-O, where the F, E, V represent the face, edge and vertex of the tetrahedron structure of CH4, while, O and 15 ACS Paragon Plus Environment


H are the oxygen and hydrogen atoms of water, respectively. The configuration opted for estimating the energies in Figure 2 is EO, which holds second highest configuration after FO. This is possibly the reason for obtaining higher energy value (5.36 kJ/mol) for the methane water interaction as compared to the average potential energies. The basis set correction factor w is calculated to have a value of 0.44, which is multiplied to BSSE and added to uncorrected interaction energies to get the CBS values. This demonstrates the validation of half-CP method ( w ≈0.50) and this proposed scheme is implemented for the interaction of water with

tetrahydrofuran and cyclopentane. For estimating the anisotropic interaction of methane with water, the potential energies are computed for angular and radial degrees of freedom. Figure 3 shows the spread of the energies as compared to the typical interaction scale, i.e. angle averaged potential. For the intermolecular distances less than 4 Å, the energies are dispersed around both the attractive and repulsive regions. As the orientations for all angular degrees of freedom in the spherical coordinates are equally probable for the guest in the water cavity, we have adopted an arithmetical average scheme for estimating the cavity potential. 1200 Average interaction energy

900 600

E CP/k (K)

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 16 of 28

300 0 -300 -600

3

4

5 6 7 8 Intermolecular distance (Å)

9

10

Figure 3. Spread of the half-counterpoise corrected interaction energies calculated at all possible angular and radial degrees of freedom. 16 ACS Paragon Plus Environment

Page 17 of 28

300

Methane 2 model fit R =0.9996 Tetrahydrofuran 2 model fit R =0.9664 Cyclopentane 2 model fit R =0.9801

200 100

E CP/k (K)

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


0

-100 -200 -300

2

3

4

5 6 7 8 9 10 Intermolecular distance (Å)

11

12

Figure 4. Angle averaged potential energy curves for interactions of various guests with water. All lines represent Kihara model fits to the arithmetically averaged potential. The possible effect of other initial orientations of guest-water for PES grid generation, as discussed in Section 2.2.1, is investigated through rotating the THF molecule by 90o along r  axis. Comparing the newly obtained cavity potential parameters (  / k = 249.37,  = 4.12) with those obtained using the original configuration (  / k = 246.75,  = 4.16), we have ignored this degree of freedom (i.e., rotation of THF molecule around r  axis) in the subsequent calculations as it does not affect the cage occupancies significantly (  L =0.9988 to 0.9991 for THF). However, for accurate hydrate phase predictions, one can consider this rotation using the average of the two sets of parameters. The angle averaged potential energies are plotted with respect to the guest-water intermolecular distances and fitted to Kihara potential curve (Figure 4). As this ab initio potentials yields negative value of hard core radius ( a ), the parameter is rather calculated using the Tee’s correlation ( 2a / (  2a )  0.1495  1.8428w ).39 The identified values of the Kihara potential function parameters for methane, THF and cyclopentane are listed in Table 1. The parameters are employed to calculate the cage occupancies of methane hydrate for both the small (  S ) and large cages (  L ) of sI hydrate and compared with the Raman spectroscopic, second 17 ACS Paragon Plus Environment


Page 18 of 28

virial coefficient and viscosity data (Table 2).37 For error analysis, an percent average absolute deviation (%AARD) method is chosen. The cage occupancies in large cages show better agreement with the experimental data with 0.33% AARD, as compared to the same predicted using the basic vdWP model with 1.33% AARD. In case of small cages, the occupancies are experimentally found to be decreasing with increase in temperature, whilst the model shows an increase. This unexpected observation is consistent with the same examined by Anderson et al.24 and Sum et al.37 The simulated ab initio cavity potential predicts the occupancies well with both the experimental data. These obtained cage occupancies are used in the vdWP model to predict the formation conditions of pure methane hydrate. Table 1 Kihara potential parameters derived from ab initio strategy Hydrate former

 (Å)

 / k (K)

a (Å)39

methane

3.6411

153.75

0.1561

tetrahydrofuran

4.1664

246.75

0.7526

cyclopentane

4.1220

250.46

0.7779

It is noteworthy that the ab initio methodology for cavity potential calculations for the promoter containing hydrates is validated first for the methane hydrate. Using the widely accepted values of reference parameters of sI hydrates, i.e.  o  1263.60, H o  - 4858.9 ,40 the hydrate formation pressures are predicted in the range of temperature (T) of 267.42–307.00 K. It is evident from Figure 5 that the model predicts the pressures close to the experimental41 values, yielding an AARD of 8.63% for 27 data points. Table 2 Experimental and predicted cage occupancies for methane hydrate

T (K)

Experimental (Raman

Experimental ( B and 

data)37

data)37

Ab initio (This work)

L

S

L

S

L

S

273.65

0.971

0.920

0.974

0.876

0.9704

0.9146

274.65

0.972

0.899

0.975

0.882

0.9718

0.9183

275.65

0.974

0.869

0.976

0.887

0.9732

0.9219

276.65

0.973

0.866

0.978

0.892

0.9746

0.9255


Page 19 of 28

160 140 Experimental Model (This work)

Pressure (MPa)

120 100 80 60 40 20 0 260

270

280 290 Temperature (K)

300

310

Figure 5. Phase equilibria prediction for methane hydrate using ab initio potential. 1.8

Raman intensity ratio Imethane/Itetrahydrofuran

6.5 6.0

1.6

5.5 1.4

5.0 4.5

1.2

4.0 1.0

3.5 3.0

nmethane/ntetrahydrofuran

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


0

1

2

3 4 5 6 Pressure (MPa)

7

8

9

0.8

Figure 6. Validating the calculated hydrate phase compositions (solid line) of CH4/THF system using ab initio potential with the Raman intensity ratio data (■).38



The in situ Raman spectroscopic analysis for the mixed sII hydrates of CH4 and THF are available by Moryama et al.24 for the relative cage occupancy quantification. In Figure 6, the experimental Raman peak ratios and the predicted hydrate phase compositions are plotted against the equilibrium pressures. It is observed that up to 3 MPa pressure, the small cage occupancies of methane hydrate steeply increase, while the variation becomes small for the higher pressures.38 It depicts that the cage occupancy ratios calculated using the proposed ab initio potential based gas hydrate model follow the similar trend of the changes in the Raman peak intensity ratios for CH4+THF binary hydrates. As the Raman intensity ratio gives a ‘qualitative’ estimation of the cage occupancy ratio of the guests in binary hydrates, the two axes in Figure 6 are chosen in such a way that the two data sets map each other. Similarly, for cyclopentane, the estimated potential parameters are used to evaluate the absolute cage occupancy of its binary hydrate with methane. The Raman spectroscopic results of Lv and Li43 suggest the occupation of methane and cyclopentane in only small and large cages, respectively. The estimation of fractional occupancies at 277.15 K results in 0.662 and 0.849 for methane, and 0.978 and 0.991 for cyclopentane at 3.5 and 12 MPa pressures, respectively. These estimated values are comparable to the reported data for cyclopentane, i.e. 0.997 and 0.995 at 3.5 and 12 MPa pressures with 277.15 K temperature. However, the estimated values of cage occupancies for methane in the binary hydrate show a significant deviation from the experimental values, i.e. 0.204 and 0.409 at 3.5 and 12 MPa pressures. This observation is consistent with that reported by Hester et al.,44 in which, they show that the Raman band areas give only the qualitative estimations of cage occupancies for the binary hydrates. This difficulty is demonstrated by Subramanian et al.42 for methane+ethane mixture hydrate. They predicted the relative hydrate composition using a calibration plot of the vapor phase compositions and ratio of the areas under the bands. However, this method has limitations for correlating effective composition ratios for the guests exhibiting in the different phases, e.g. CH4 in vapor and THF in liquid. Noting this, we have adopted the method of Moryama et al.24 to correlate Raman intensities directly to the molar ratios for relative estimation of cage occupancies.24 A recent X-ray diffraction study on lattice distortion caused by THF and cyclopentane reported their cell sizes of 17.1740 and 17.2464, respectively, even though their molecular sizes are similar.45 This is due to the different lattice stabilization energies for these two promotors. 20 ACS Paragon Plus Environment

Page 20 of 28

Page 21 of 28 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


However, this does not affect the ab initio potentials obtained in this study as they are a function of intermolecular distance between the guest and water. These illustrations show that that the strategy developed for the promoter containing gas hydrates is equally applicable to the small guests, e.g. methane, and thus providing a generalized ground to estimate cavity potential. During the simulation, the possible computational deficiencies can be handled by forming denser grid for potential energy surface; selecting sophisticated ab initio methods and opting proper potential model for fitting to ab initio energies. Table 3 Comparison of QM methods and basis sets applied to estimate cavity potential Reference

Processor used

QM method /Basis set Time#

Time$

Guests

Cao et al.19-21

180 MHz,

MP2/6-31++G(2d,2p)

min

4.1 sec*

CH4

Silicon

MP2/cc-pVQZ

days

4 min

CH4

MP2/6-31++G(3d,3p)

~1 day

8.2 sec*

C1-C3, CO2

graphics Klauda and

195 MHz

Sandler22,23

MIPS R10000

Anderson et

N.A.

MP2/aug-cc-pvQZ

N.A.

22.5 min

CH4, Ar

N.A.

MP2

N.A.

N.A.

CH4, CO2

N.A.

MP2/aug-cc-pvTZ

~10 min

1.3 min

CO2

2.5 GHz AMD

SCS-MP2/aug-cc-

Computing Anisotropic Cavity Potential for Clathrate Hydrates - The

Recommend Documents