CC BY
71
УДК 548.562
Г.С. Смирнов, В.В. Стегайлов
Молекулярно-динамические модели газовых гидратов и описание фазовых диаграмм
Ключевые слова:
газовые гидраты,
молекулярное
моделирование,
фазовые
диаграммы,
межмолекулярные
взаимодействия.
Keywords:
gas hydrates, molecular modelling, phase diagrams, intermolecular interactions.
Гидраты природных газов или газовые гидраты - это нестехиометрические соединения включения, в которых молекулы газа (молекулы-гости) заключены в полостях трехмерной решетки из молекул воды (каркас хозяина). Стабилизация водных кла-тратных каркасов, термодинамически менее стабильных, чем лед или жидкая вода при тех же условиях, обеспечивается за счет ван-дер-ваальсовых взаимодействий «гость - хозяин» [1].
Месторождения природных газовых гидратов встречаются на шельфе и в зоне вечной мерзлоты. Они активно изучаются уже более сорока лет, так как содержат большое количество метана и рассматриваются как перспективный источник топлива. Вместе с тем газовые гидраты, существующие при температурах и давлениях, близких к условиям их разложения, представляют экологическую опасность - в случае смещения теплового равновесия освободившийся метан может внести значительный негативный вклад в парниковый эффект. Поэтому знание динамических, термодинамических и механических свойств газовых гидратов важно для оценки их поведения при разработке газогидратных месторождений.
За последние годы было проведено множество теоретических, экспериментальных и компьютерных исследований газовых гидратов. Так, исследования месторождений на Аляске, в Канаде и Японии показывают, что добыча метана из природных гидратов может быть выполнена с использованием существующих технологий. Более того, около побережья Японии недавно была произведена первая в мире пробная добыча газа из месторождения в желобе Нанкай [2].
Молекулярная динамика позволяет исследовать свойства вещества на микроскопическом уровне, поэтому ее развитие привело к пониманию ряда фундаментальных свойств газовых гидратов. Были исследованы как механические, так и термодинамические свойства газовых гидратов в широком диапазоне температур и давлений для различных газовых гидратов [3-9], достигнуто определенное понимание процессов гидратообразования [10, 11].
Хотя кристаллы газовых гидратов по большей части состоят из молекул воды, их структура и свойства существенно отличаются от любых известных льдов. В отличие ото льдов, которые формируются из чистого вещества, газовые гидраты не могут образовываться без молекул газа определенного размера. В природе наиболее распространены три структуры газовых гидратов - простая кубическая (КС-1), гранецентри-рованная кубическая (КС-11) и гексагональная (ГС-Ш). В «классических» клатратных гидратах молекулы воды формируют полости в форме многогранников, в которых находятся молекулы газа. Для их описания применяется обозначение пт, где п - число ребер в грани, а т - число граней с п ребрами. Например, в элементарной ячейке структуры КС-1 содержатся 46 молекул воды, которые образуют два двенадцатигранника 512 и шесть четырнадцатигранников 51262.
В области высоких давлений молекулы газа небольшого размера могут образовывать соединения включения с различными формами водного льда. Данные структуры представляют интерес с фундаментальной точки зрения, а также важны для понимания строения ледяных планет и лун.
Полости в газовых гидратах могут быть заполнены как полностью, так и частично. В области высоких давлений также возможно включение в полости двух и более молекул. Зная размеры полостей и гостевых молекул, можно на качественном уровне предсказывать возможность образования тех или иных клатрат-ных структур. Так, КС-1 образуется молекулами размерами от 4,2 до 6 А (примером могут служить структуры СН4, С2Н6, С02, И28). Маленькие (размерами меньше 4,2 А) и большие (6 < й < 7 А) молекулы формируют КС-11. Еще большие молекулы (7 < й < 9 А), например изопентан или неогексан, могут формировать гексагональную структуру, заполняя большие полости. Для образования термодинамически устойчивых соединений при этом необходимо заполнение малых полостей. Они заполняются другими молекулами, например метаном, сульфидом водорода или азотом. Приведенные рассуждения справедливы только при давлениях порядка 1-100 МПа. При более высоких давлениях возможна реализация нескольких газогид-ратных структур в одной системе. Например, в системах «аргон - вода» и «водород - вода» обнаружено не менее четырех типов гидратов.
Для разработки месторождений гидратов необходимо знание их фазовой диаграммы. Существующие модели хорошо описывают диаграмму чистого гидрата метана в широком диапазоне давлений и температур. Однако на границу устойчивости гидратов могут влиять параметры окружающей среды, роль которых изучена хуже. Например, в пористой среде с диаметром пор меньше 55 нм наблюдается сдвиг кривой фазовой равновесия [12]. Аналогично на кривую фазового равновесия влияет наличие примесных солей в воде, что критично при разработке морских месторождений гидратов. Указанные эффекты могут быть учтены при построении многомасштабных моделей формирования и разрушения гидратов. Одним из этапов такого моделирования может выступать молекулярно-динамическое моделирование кривой плавления газовых гидратов, рассматриваемое в данной работе.
Молекулярно-динамическая модель
В основе метода молекулярной динамики лежит численное интегрирование классических уравнений движения Ньютона для системы N частиц. Свойства системы закладываются в потенциал межатомного взаимодействия. Для
расчета взаимодействия между молекулами были выбраны модели, основанные на потенциалах Леннарда-Джонса и Кулона:
и=ч
г..
v ч у
ч *
г
V ч
Коэффициенты е, с и q в потенциале выбираются для согласования с экспериментальными данными. В данной работе используются следующие параметризации: для воды - 8РС/Е [13], Т1Р4Р/2005 [14] и Т1Р4Р/1СЕ [15]; для метана - из работ [16, 17]. Перекрестное взаимодействие задавалось правилами Лоренца-Бертло:
Ст.. + ст..
СТЧ=-
=хут^ч! .
Вышеуказанные потенциалы были разработаны для чистых веществ, в то же время эмпирические правила Лоренца-Бертло позволяют применять потенциалы и для газовых гидратов. Варьирование параметра х позволяет добиться лучшего согласия с экспериментом. В статье [18] рассмотрен разработанный потенциал с х = 1,07 для модели воды Т1Р4Р/2005, позволивший неявным образом учесть поляризационные эффекты взаимодействия молекул метана с клатратной структурой.
Расчет кривой плавления КС-! гидрата метана
Для расчета кривой плавления гидрата метана рассматривался параллелепипед, в одной половине которого находился кристалл, а в другой - смесь жидкой воды и метана. Размеры системы составляли от 2*2*4 до 6*6*12 элементарных ячеек (обсуждение размерного эффекта приведено ниже). Далее в ходе короткого NPT-расчета1 система выводилась в желаемый диапазон температур и давлений. Для определения равновесных значений давления и температуры проводился молекулярно-динамический расчет в NVE-ансамбле2, в ходе которого система выходила на равновесие [19].
В зависимости от начальной кинетической энергии возможны два процесса - рост
1 Моделирование системы с фиксированными количеством частиц давлением (Р)
и температурой (Т).
2 Моделирование системы с постоянными значениями числа частиц объема (V) и энергии (Е).
кристалла либо его плавление. Если начальная температура системы слишком высока (или, наоборот, низка), кристалл полностью расплавится (или, наоборот, вся жидкость кристаллизуется). Такой процесс не позволит определить точку фазового равновесия, поскольку в процессе расчета должна образоваться равновесная трехфазная система «кристалл - жидкость - газ». Именно поэтому степени перегрева или переохлаждения не должны быть слишком большими. Пример КС-1 гидрата метана показан на рис. 1.
Анализ показывает, что молекулы метана собираются в пузырь цилиндрической формы. Это затрудняет рост кристалла, а поверхностное натяжение может влиять на положение кривой фазового равновесия. Для уменьшения времени расчета исследовались только перегретые состояния, так как кинетика процесса плавления быстрее процесса кристаллизации.
Аналогичные расчеты кривой фазового равновесия приведены в работе М. Конде и К. Веги [20]. Они проводили длительные (до 1 мкс) вычисления молекулярно-динамических траекторий в ЛУТ-ансамбле. В начале расчета в одной половине прямоугольной расчетной ячейки находится кристалл гидрата метана, в другой - его расплав. Поскольку точка фазового равновесия заранее неизвестна, в системе наблюдается рост кристалла либо его плавление. Такой процесс зависит от кинетической энергии всех молекул.
Направления водородных связей между молекулами воды в клатратном каркасе (сами молекулы воды не показаны)
Рис. 1. Визуализация молекулярной модели КС-1 гидрата метана
Поскольку полная энергия не является постоянной величиной, потенциальная энергия изменяется при фазовом переходе. Уменьшение потенциальной энергии соответствует росту кристалла, а ее рост - плавлению. Проводя серию расчетов с различной температурой, но постоянным давлением, можно определить равновесную температуру как среднее между наибольшей Т, при которой наблюдается рост кристалла, и наименьшей Т, при которой наблюдается плавление. Главным недостатком этого метода является необходимость проведения большого количества длительных расчетов, что затрудняет проверку зависимости результата от числа частиц.
В ходе расчетов выяснилось, что с увеличением размера системы увеличивается значение температуры плавления. Подобный эффект может быть вызван наличием поверхности раздела фаз, свободная энергия которой может сдвигать точку плавления. Разница становится несущественной при использовании пяти элементарных ячеек по каждому направлению в поперечном сечении, именно такая система исследовалась в дальнейших расчетах. Результаты расчетов для разных моделей воды приведены на рис. 2 [19].
Все потенциалы качественно воспроизводят вид кривой плавления, хотя в области высоких давлений наблюдается расхождение с экспериментом. Тем не менее, приведенные примеры расчетов кривой плавления КС-1 гидрата метана показывают, что кривые фазовых равновесий газовых гидратов могут изучаться с использованием методов классической молекулярной динамики с достаточно высокой степенью точности (погрешность по температуре составляет ~ 20 К).
Оценка точности описания фазовых диаграмм
Атомистические модели для газовых гидратов существенным образом основаны на моделях межатомных и межмолекулярных взаимодействий. Однако описание нековалентных взаимодействий из первых принципов представляет собой трудную задачу. Квантово-механические «пост-хартри-фоковские» методы сегодня только показывают будущие перспективы построения фазовых диаграмм на молекулярном уровне без учета эмпирических поправок для дисперсионного взаимодействия [23]. Хорошо развитые методы теории функци-
<2 800
1600
400
200
\ \
А □ \ <> 1 1 1 1 ( > I
— экспериментальнаялинияплавления [21] ФН^ данные модельных расчетов работы [20] данные модельных расчетов работы [22] ОПО^ результаты авторов данной статьи для тех же моделей межмолекулярного взаимодействия
/ - 1 / / / / / -/—
А □ - - " ----- . - - "' у о
260
280
300
320 340
Температура, К
Рис. 2. Фазовая диаграмма КС-1 гидрата метана:
синие символы - модель воды Т1Р4Р/1се; красные - Т1Р4Р/2005 (х = 1); зеленые - Т1Р4Р/2005 (х = 1,07); черные - 8РС/Б; область между пунктирными линиями (интерполяциями результатов расчетов в рамках моделей Т1Р4Р/1се и Т1Р4Р/2005 (X = 1,07)) иллюстрирует текущую точность определения кривой плавления по результатам моделирования
0
онала электронной плотности не могут описывать дисперсионное взаимодействие, но учитывают эффекты перераспределения электронной плотности и поляризации для моделей ограниченного размера [24, 25]. Поэтому на сегодняшний день квантовые расчеты электронной структуры и первопринципная молекулярная динамика для супрамолекулярных систем не имеют существенных преимуществ по сравнению с классической молекулярной динамикой, основанной на эмпирических потенциалах взаимодействия.
Молекулярно-динамические модели для клатратных гидратов существенным образом зависят от точности описания взаимодействия между молекулами воды. Точность таких моделей сильно возросла за последнее десятилетие. Наиболее сложные модели, учитывающие поляризацию, лучше описывают структурные свойства [26] и растворимость молекул газа [27], тогда как модели без поляризации лучше подходят для описания фазовой диаграммы [28, 29].
Для оценки систематических ошибок при описании фазовых диаграмм можно применять эвристический подход, основанный на законе соответствующих состояний, который хорошо работает для простых металлов, жидкостей и газов. Взяв пять экспериментальных и
рассчитанных тройных точек на фазовой диаграмме льда (рис. 3), находим аффинное преобразование А, наилучшим образом переводящее
400
+ VI + + ♦ :
+ + V • ' * * . ■ + + ♦ •
♦ + тт ♦ ♦♦♦♦♦«■ I V-/ .V ♦♦♦А: \
•
200
240
280
Температура, К
экспериментальная фазовая диаграмма воды
(1-У - типы льдов)
расчетная фазовая диаграмма
для Т1Р4Р/1се модели
фазовая диаграмма после аффинного
преобразования^
пять тройных точек
Рис. 3. Фазовая диаграмма, построенная после аффинного преобразования А, полученного по пяти тройным точкам
а 800
0
рассчитанную фазовую диаграмму в экспериментальную [30]:
( pTIPAP/Ice \
TT
1
= A
' pExp Tt Exp
1
V У
Используя обратное преобразование, можно найти систематическую ошибку описания фазовых диаграмм газогидратных структур. Анализ показывает, что систематическая погрешность составляет порядка ±0,2 ГПа и ±50 К [30].
В работе приведен краткий обзор свойств и приложений газовых гидратов и исследования их фазовых диаграмм методами молекулярно-
динамического моделирования. Граница трехфазного равновесия КС-1 гидрата метана определена для различных моделей воды. Установлено, что модель воды ТГР4РЛсе дает завышенное значение границы трехфазного равновесия, остальные модели - заниженное. Точность описания кривой плавления составляет 20 К. Проанализирована точность современных атомистических моделей межмолекулярных потенциалов для воды как основы описания газовых гидратов. Показано, что порядок систематической погрешности составляет десятые доли ГПа по давлению и десятки градусов по температуре.
Работа выполнена при частичной поддержке по гранту РФФИ 14-08-31550 мол а.
Список литературы
1. Sloan E.D. Clathrate hydrates of natural gases / E.D. Sloan, C. Koh. - 3ed ed. - CRC Press, 2007. -701 c.
2. Yamamoto K. Operational overview of the first offshore production test of methane hydrates in the Eastern Nankai Trough / K. Yamamoto,
Y. Terao, T. Fujii et al. // Offshore technol. conf. -2014. - P. 1-11.
3. Westacott R. Full-coordinate free-energy minimisation for complex molecular crystals: type I hydrates / R. Westacott, P. Rodger // Chem. Phys. Lett. - 1996. - V. 262. - № 1-2. - P. 47-51.
4. Moon C. Molecular dynamics study of gas hydrate formation / C. Moon, P.C. Taylor, P.M. Rodger // J. Am. Chem. Soc. - 2003. - V. 125. - № 16. -
P. 4706-7.
5. English N.J. Theoretical studies of the kinetics of methane hydrate crystallization
in external electromagnetic fields / N.J. English, J.M.D. MacElroy // J. Chem. Phys. - 2004. -V. 120. - № 21. - P. 10247-56.
6. English N.J. Molecular-dynamics simulations of methane hydrate dissociation / N.J. English, J.K. Johnson, C.E. Taylor // J. Chem. Phys. -2005. - V. 123. - № 24. - P. 244503.
7. Istomin V.A. Metastable states of gas hydrates / V.A. Istomin, V.G. Kvon, V.A. Durov // Gas Ind. Russ. - 2006. - V. 4. - P. 13-16.
8. Tung Y.-T. The growth of structure I methane hydrate from molecular dynamics simulations / Y.-T. Tung, L.-J. Chen, Y.-P. Chen, S.-T. Lin // J. Phys. Chem. B. - 2010. - V. 114. - № 33. -P. 10804-13.
9. Walsh M.R. Methane hydrate nucleation rates from molecular dynamics simulations: effects of aqueous methane concentration, interfacial curvature, and system size / M.R. Walsh,
G.T. Beckham, C. Koh et al. // J. Phys. Chem. C. -2011. - V. 115. - № 43. - P. 21241-21248.
10. Jacobson L.C. Nucleation pathways of clathrate hydrates: effect of guest size and solubility / L.C. Jacobson, W. Hujo, V. Molinero // J. Phys. Chem. B. - 2010. - V. 114. - № 43 - P. 13796807.
11. Knott B.C. Homogeneous nucleation of methane hydrates: unrealistic under realistic conditions / B.C. Knott, V. Molinero, M.F. Doherty et al. //
J. Am. Chem. Soc. - 2012. - V. 134. - № 48. -P. 19544-7.
12. Turner D.J. Sensitivity of methane hydrate phase equilibria to sediment pore size / D.J. Turner, R.S. Cherry, E.D. Sloan // Fluid Phase Equilib. -2005. - V. 228-229. - P. 505-510.
13. Berendsen H.J.C. The missing term in effective pair potentials / H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma // J. Phys. Chem. - 1987. - V. 91 -№ 24. - P. 6269-6271.
14. Abascal J.L.F. A general purpose model
for the condensed phases of water: TIP4P/2005 / J.L.F. Abascal, C. Vega // J. Chem. Phys. - 2005. -V. 123. - № 23. - P. 234505.
15. Abascal J.L.F. A potential model for the study of ices and amorphous water: TIP4P/Ice / J.L.F. Abascal, E. Sanz, R. García Fernández
et al. // J. Chem. Phys. - 2005. - V. 122. - № 23. -P. 234511.
16. Guillot B. A computer simulation study
of the temperature dependence of the hydrophobic hydration / B. Guillot, Y. Guissani // J. Chem. Phys. - 1993. - V. 99. - № 10. - P. 8075.
17. Paschek D. Temperature dependence
of the hydrophobic hydration and interaction of simple solutes: an examination of five popular water models / D. Paschek // J. Chem. Phys. -2004. - V. 120. - № 14. - P. 6674-90.
18. Docherty H. A potential model for methane in water describing correctly the solubility of the gas and the properties of the methane hydrate / H. Docherty, A. Galindo, C. Vega
et al. // J. Chem. Phys. - 2006. - V. 125. - № 7. -P. 074510.
19. Smirnov G.S. Melting and superheating
of si methane hydrate: molecular dynamics study / G.S. Smirnov, V.V. Stegailov // J. Chem. Phys. -2012. - V. 136. - № 4. - P. 044523.
20. Conde M.M. Determining the three-phase coexistence line in methane hydrates using computer simulations / M.M. Conde, C. Vega // J. Chem. Phys. - 2010. - V. 133. - № 6. -
P. 064507.
21. Dyadin Y.A. Decomposition of methane hydrates up to 15 kbar / Y.A. Dyadin, E.Y. Aladko,
E.G. Larionov // Mendeleev Commun. - 1997. -V. 7 - № 1. - P. 34-35.
22. Jensen L. Calculation of liquid water-hydrate-methane vapor phase equilibria from molecular simulations / L. Jensen, K. Thomsen, N. von Solms et al. // J. Phys. Chem. B. - 2010. -V. 114. - № 17. - P. 5775-82.
23. Ben M.Del. Second-order Maller-Plesset perturbation theory in the condensed phase: An efficient and massively parallel gaussian
and plane waves approach / M.Del. Ben, J. Hutter, J. Vandevondele // J. Chem. Theory Comput. -2012. - V. 8. - № 11. - P. 4177-4188.
24. Wang J. Molecular-dynamics and first-principles calculations of raman spectra and molecular and electronic structure of hydrogen clusters
in hydrogen clathrate hydrate / J. Wang, H. Lu, J.A. Ripmeester et al. // J. Phys. Chem. C. -2010. - V. 114. - № 49 - P. 21042-21050.
25. Zhang J. First principles molecular dynamics study of filled ice hydrogen hydrate / J. Zhang, J.-L. Kuo, T. Iitaka // J. Chem. Phys. - 2012. -V. 137. - № 8. - P. 084505.
26. Jiang H. Molecular dynamics simulations
of methane hydrate using polarizable force fields / H. Jiang, K.D. Jordan, C.E. Taylor // J. Phys. Chem. B. - 2007. - V. 111. - № 23. - P. 6486-92.
27. Dyer P.J. The importance of polarizability
in the modeling of solubility: quantifying the effect of solute polarizability on the solubility of small nonpolar solutes in popular models of water / P. J. Dyer, H. Docherty, P.T. Cummings // J. Chem. Phys. - 2008. - V. 129. - № 2. - P. 024508.
28. Kiss P.T. Testing recent charge-on-spring type polarizable water models. I. Melting temperature and ice properties / P.T. Kiss, P. Bertsyk,
A. Baranyai // J. Chem. Phys. - 2012. - V. 137. -№ 19. - P. 194102.
29. Conde M.M. Note: A simple correlation to locate the three phase coexistence line in methanehydrate simulations / M.M. Conde, C. Vega //
J. Chem. Phys. - 2013. - V. 138. - № 5. -P. 056101.
30. Smirnov G.S. Toward determination of the new hydrogen hydrate clathrate structures /
G.S. Smirnov, V.V. Stegailov // J. Phys. Chem. Lett. - 2013. - V. 4. - № 21. - P. 3560-3564.
