Study for the stability of CO2 clathrate-hydrate using molecular dynamics simulation

Energy Convers. Mgmt Vol. 37, Nos 6-8, pp. 1087-1092, 1996

Pergamon

0196-8904(95)00301-0

Copyright 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0196-8904/96 $15.00 + 0.00

S T U D Y F O R T H E STABILITY OF CO2 C L A T H R A T E - H Y D R A T E MOLECULAR

USING

DYNAMICS SIMULATION

S. HIRAI t, K. OKAZAKI, S. KURAOKA and K. KAWAMURA tt Research Center for Carbon Recycling and Utilization, ~*Department of Earth and Planetary Sciences, Faculty of Science, Tokyo Institute of Technology, Ohokayama, Meguroku, Tokyo, 152,

Japan

Abstract - Investigation using molecular dynamics simulation has been conducted in order to clarify the stability and to obtain the fundamental understanding of CO2 clathrate-hydrate. It was revealed from the MD simulation that the CO2 clathrate-hydrate is unstable as compared with both empty and Argon clathratehydrates. The reason for the unstableness is discussed based on the database obtained from the simulation. The repulsive force acting between the O atoms of CO2 and O atom in H20 consisting the cage have a destabilizing effect on the CO2 clathrate-hydrate lattice structure. 1. INTRODUCTION Clathrate-hydrate is a crystalline solid where water molecules are linked through hydrogen bonding making lattice structure with cavities in which inclusion of second kind molecules (guest molecules) in the cavities is observed. When the guest molecule is carbon dioxide (CO2), it is called a CO2 clathrate-hydrate. Two types of clathrates are found, type I and II ,and CO2 clathrate-hydrate belongs to type I where type I consists of large tetrakaidechedral cages and small cages. CO2 molecules are only contained in the large cages due to its large molecular size. CO2 clathrate-hydrate is formed under the conditions of temperature less than 10 °C and pressure higher than 44 atm. Recently, CO2 clathrate-hydrate has been attracting special interest. CO2 sequestration in ocean, which are (1) dissolution of liquefied CO2 at intermediate depth of sea and (2) disposal of liquid CO2 at seafloor, is considered to be one of mitigation strategy in order to avoid the global warming due to increase of CO2 concentration in the atmosphere [1,2]. CO2 clathrate-hydrate film is formed at the interface between sea water and liquid CO2 and CO2 clathrate-hydrate film is considered to have a large effect on the diffusion of CO2 into the seawater. The understanding of the CO2 clathrate-hydrate, particularly its stability, is remarkably important, because CO2 clathrate-hydrate film acts as a resistance layer for the dissolution of CO2. Besides the experimental study of CO2 clathrate-hydrate [3,4], molecular dynamics simulation study provides a numerous data to consider a relation between its characteristic lattice structure and stability from molecular viewpoints. Tse et al made computer simulation studies of the clathrate-hydrate of methane [5] and xenon [6]. Marchi and Mountain investigated the stability of clathrate crystal containing krypton relative to empty hydrate [7]. In author's knowledge, however, there has been no reports for the molecular dynamics simulation study of CO2 clathrate-hydrate. The object of the present paper is to investigate the stability of the CO2 clathrate-hydrate using the database obtained from the molecular dynamics simulation. It was made clear that CO2 clathrate-hydrate is unstable as compared with both empty and Ar clathrate-hydrates. The interaction between the encaged tTo whom all correspondence should be addressed. 1087

1088

HIRAI et al.: STABILITY OF CO2 CLATHRATE-HYDRATE

guest CO2 molecules and the water molecules making a lattice structure was elucidated. 2. THE INTERATOMIC POTENTIAL MODEL AND THE SIMULATION METHOD Many rigid molecular models have been made for molecular simulation of H20. Though these models possess ability to reproduce the characteristics of liquid water, some difficulties are considered to appear in the simulation of solid ice. The present calculation applies a elaborate model [8] which treats Hz0 systems as an assemblage of hydrogen and oxygen atoms. The model expresses the inter- and intramolecular interactions in terms of a single interatomic potential functions and motions of H and 0 atoms have a total freedom. The interatomic potential function is made of two body and three body terms. The two body terms are expressed by eq. (1) for pairs of atoms [8] .

ui,

=

‘izie2 -+f,(& 'ij

1

‘“ids,“’ -c,cI rti6

+b,)exp i

1

1

+f,D,(ev[-W,,(r,

-r~*)]-2eXP[-Aj(ri,

(1)

-rij*)])

Here, ‘ij is a interatomic distance, f, a constant for unit adaptations for these terms. The first, the second and the third terms represents Coulomb, short range repulsions and van der Waals interaction terms, respectively. The fourth term (the Morse function) are only employed to estimate the interactions between O-H atoms and C-O atoms. The three body part is only employed for H-O-H groups as LiHoH(e”OII)=-2~~(cos[2(eHOH -%)]-l]@-

(2)

where oHoHis the angle of H-O-H and k, and k, define the effective range of three body potential by the following equation k, =I/{exp[g,(r,,(i)-r,,,)]+l)

(3)

The three body forces apply perpendicular to the O-H bond and its magnitude is given by PH(0 = -2f, sin[2(&,, - &)]~/ro,W

(4)

The parameters that appear in above equations are listed in Table 1. The time advance scheme applied in the present calculation is that of Verlet algorithm and Ewald Table 1. Parameters of interatomic potential models

. 1.841:0.036:

1.885:0.690:

1.86

HIRAI et al.: STABILITYOF CO2 CLATHRATE-HYDRATE

1089

summation is applied to evaluate Coulomb forces. The basic cell is a cubic one where 6 large cages and 2 small cages are involved. The time integration was carded out using At = 0.4 fs under the constant temperature and pressure conditions. The initial coordinates of the atoms are determined by the crystal structure obtained by the neutron-ray diffraction experiments [9]. 3. RESULTS AND DISCUSSIONS Figs. l(a) and (b) illustrate the profile of M.S.D.(mean-square displacement) profiles of O atoms in I-I20 consisting the cages for CO2 clathrate-hydrate and empty one (no guest molecule in the cage), respectively. Note that the M.S.D. is estimated using the O atoms in H20, which are those making the cage structure. M.S.D. is evaluated from the following equation.

M.S.D.= ~, {R,(t)- R,(to)}2/N

(5)

Here, Ri are the coordinates of O atoms consisting the cages (O atoms in H20 molecules) and N is the total number of the O atoms. M.S.D. of CO2 clathrate-hydrate (Fig. 1(a)) shows large value from the beginning of the calculation whereas empty clathrate-hydrate (Fig. l(b)) slowly reaches its equilibrium state. In addition, the magnitude of the M.S.D. of CO2 clathrate-hydrate in the equilibrium state is larger than that of empty one. These are considered to be caused by the effect of the CO2 molecules contained in the large cage, where initial condition of the calculation is obtained by the crystalline structure of clathrate-hydrate without CO2 molecules. Figs.1 (a) and (b) also show that the amplitude of M.S.D. of CO2 clathrate-hydrate in the equilibrium states is larger than that of empty one. The present calculation controls the lattice parameters for constant pressure condition and its effect is included in M.S.D. Figs. 1 (a) and (b) compare the clathrate-hydrate between CO2 and empty ones in the same calculation conditions, where it is considered that the qualitative estimation for the effect of CO2 molecule could be made. The trajectory of a CO2 molecule inside a large cage is shown in Fig. 2(a) for the temperature 250K. The relation between the cage structure and the coordinates are shown in Fig. 2(b). The CO2 molecule shows complicated three dimensional behavior but it should be noted that the position of C atom in CO2 is restricted to the center part of the large cage and two O atoms of CO2 rotates around the C atom. In order to investigate the stability of the CO2 clathrate-hydrate, we focus on the pair-correlation function (P.C.F.(r)) between O-O molecules consisting the cage (not the O molecules of CO2). P.C.F.(r) can be evaluated from the following equation. P .C . F. ( .r ) .= n / [

I

I

I

4]r(N, , ,NI/VIr,2Ar .•,

]

(6)

i

*~< a"0'2I''' 05

~< 0.2 r~

o5

0.1

0.1

1'o

2'o

Timeps

3'0

(a) CO2 clathrate-hydrate

o

¢o

2'o

3'o

Timeps (b) Emptyclathrate-hydrate

Fig. 1 M.S.D. profiles of O atoms in H20 consisting the cage ( Temperature 250 K )

HIRAI et al.: STABILITY OF CO2 CLATHRATE-HYDRATE

1090 X

11

12

13

z

10

(a) Temperature 250 K

(b) The relation between the cage structure and the coordinates

Fig. 2. Trajectory of a CO2 molecule inside a large cage (

4-,

I ”

1 ‘.

”

”

1 * ”

4

0 :C atom, 0 I

I

----. Ar L6 . --_-----_ Empty CO2 L6-

I

I.

3.

2.5 Distance A (a> 1st. peak

I c..

3

z

15

and 0 : 0 atoms )

’

,

.

----Ar L6 -_------- Empty CO2 L6-

I

4

I

I

5

Distance A

(b) 2nd. peak

Fig. 3. P.C.F. profiles of CO??,empty and Ar clathrate-hydrate (Temperature 250K) Here, nij is the number of pairs of 0 atoms that exists in the distance from r, -Ar/2 to r, +Ar/2. Ni and Nj are the number of 0 atoms of Hz0 in the basic cell and V is the volume of the basic cell. Higher peak distribution of P.C.F.(r) represents the structure is stable and broad one shows unstable. A comparison of the first and the second peak profiles of P.C.F.(r) of CO2 clathrate-hydrate with the empty and Ar clathrate-hydrate for temperature 250K and 1OOKare shown in Figs. 3. It can be seen from the Figs. 3 that the CO2 clathrate-hydrate is unstable. It should be noted that Ar encaged in the cage serve to stabilize the cage structure as compared with that of the empty clathrate-hydrate. Therefore, guest molecules enclosed in the cages has counter effects, stabilizing effect in the case of Ar and destabilizing effect in the case of CO2, as compared with empty clathrate. Thus, CO2 molecule has a large destabilizing effect on the structure of clathrate-hydrate. We now consider its mechanism in detail. In Fig. 4, time variation of the distance between 0 atoms (marked with arrows in Figs. 5) which are placed mutually at the opposite side of the large cage is shown. The O-O distance shows the stretching vibration behavior of large cage in the O-O direction. The maximum and minimum values of O-O distance are around 9.4A and 7.8A, respectively. Figs. 5 (a), (b) illustrate relation between the direction of the axis of CO2 molecule and the large cage structure in which CO2 molecule is involved when the O-O distance is in the maximum value (9.4, A& in Fig.4 ) and minimum value (7.8, A ‘l’in Fig.4 ), respectively. It can be recognized from the figure that when the cage is extended (O-O distance

HIRAI et

STABILITYOF CO2 CLATHRATE-HYDRATE

al.: .

.

.

.

1091

I

.< 9

t5

o,

© 8

0

1000 2000 Time fs Fig. 4. Time variation of the distance between O atoms placed mutually opposite side of the cage ( Temperature 250 K ) 9.4A), the O-O axis is nearly parallel to the CO2 axis whereas when the cage is contracted (O-O distance 7.8/~), the O-O axis is nearly perpendicular to the CO2 axis. Thus, stretching vibration behavior of large cage shows strong correlation to the axis of CO2 molecule contained. CO2 molecule is making a three dimensional motion in the large cage and Fig.6 shows the time variation of the minimum distance between O atom in CO2 and O atom consisting the large cage. Its value shows variation between 3A - 5/~ and its value can be correlated with the potential energy profile between these O atoms in Fig.7. U is the summation of Coulomb, short range repulsions and van der Waals potential energy. Fig. 7 shows that, in the range of 3/~ - 5]k, Coulomb force is the dominating factor of total potential energy and it acts as a repulsion effect between the O atoms ( O atom in CO2 and O atom consisting the large cage). The clathrate cage also include H atoms ( in H20 molecules ) and Coulomb force between H atoms in H20 and O atom in CO2 acts as an attractive force. Due to the contradictory effects between the O atom and the H atoms in H20 by the O atom in CO2, the hydrogen bondings between the H20 molecules are also weakened and the clathrate structure is also destabilized.

(a) O-O (marked with arrow) distance is max.

(b) O-O (marked with arrow) distance is min.

Fig. 5. Relation between the axis of a CO2 molecule and cage structure (Temperature 250 K, White balls : O atoms in H20, gray balls : CO2)

1092

HIRAI et al.: STABILITY OF CO2 C L A T H R A T E - H Y D R A T E

'•

[xl 0 -lr] 1

U1 (Coulomb) U2(Short range repulsions) U3(van der Waals) Usum.(Ul+U2+U3)

...... t~ e--

4

(3)

O

Q.

0 c-

[/I

a i

i

I

1000

t

=

i

i

2000

T i m e fs

Fig. 6. Time variation of minimum distance between O atom(in CO2) and O atom (in H20)

'

~

'

~

'

4

~

6

Distance A Fig. 7. Potential energy profile of O atom (in CO2) andO atom (in H20)

4. CONCLUSION Molecular dynamics simulation study of CO2 clathrate-hydrate has been conducted to investigate the relation between its stability and characteristic lattice structure. It was revealed from the biD simulation that the CO2 clathrate-hydrate is unstable as compared with both empty and Argon clathrate-hydrates. The behavior of C atom in CO2 is restricted to the center part of the cage and O atoms in CO2 rotates around the C atom. The unstableness of the cage is caused by the interaction between O atom in CO2 and O atom consisting the large cage. In addition, the Coulomb force among these O atoms is the dominating factor that destabilizes the CO2 clathrate-hydrate lattice structure. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9.

C. Marchetti, Climate Change, 1, 59 (1977). M. Steinberg, et al., Brookhaven National Laboratory Report, OE/CH/IX)016 (1984). I. Aya, K. Yamane, and N. Yamada, Trans. Japan Society Mech. Eng., 59, 1210 (1993). S. Hirai, K. Okazaki, N. Araki, K. Yoshimoto, H. Ito, and K. Hijikata, Int. Conf. Carbon Dioxide Removal 2, (submitted to Energy Convers. Mgmt.) J.S. Tse, M. L. Klein, and I. R. McDonald, J. Phys. Chem., 87, 4198 (1983). J.S. Tse, M. L. Klein, and I. R. McDonald, J. Chem. Phys., 78, 2096 (1983). M. Marchi, and R. D. Mountain, J. Chem. Phys. ,86, 6454 (1987). N. Kumagai, K. Kawarnura and T. Yokokawa, Molecular Simulation., 12, 177 (1994). F. Hollander and G. A. Jeffrey, J. Chem. Phys., 66,4699 (1977).