CC BY
59
УДК 541.182.43:532.65
Вестник СЦбГУ. Сер. 4, 2004, вып. 4
Г. В. Мрджикова, Е. Н. Бродская
МОДЕЛИРОВАНИЕ МАЛЫХ КЛАСТЕРОВ ВОДЫ В НЕПОЛЯРНОЙ СРЕДЕ Н-ГЕКСАНА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
Высокодисперсные эмульсии масло-вода занимают особое место в природе, что обусловливает их большое практическое и теоретическое значение. Межфазная поверхность между двумя несмешивающимися жидкостями представляет собой трудно исследуемый объект. Развитие экспериментальных методов и большие вычислительные возможности современных компьютеров позволили разработать широкий спектр методик, используемых для анализа таких систем. При этом очень часто обращаются к методам численного моделирования, хорошо зарекомендовавшим себя в рассмотрении деталей структуры и динамики поверхностных слоев.
В одной из первых работ, посвященных исследованию систем вода-слабо полярная жидкость [1], методом молекулярной динамики (МД) моделировали поверхностный слой между водой й дихлорэтаном при 300 К. При оценке толщины поверхностного слоя получили значения в несколько ангстрем. Примерно в двух-трех молекулярных слоях наблюдали преимущественную ориентацию дипольных моментов молекул воды параллельно поверхности. На основе формулы Кирквуда-Баффа было рассчитано поверхностное натяжение, равное 35 ±20 мНм"1, достаточно хорошо согласующееся с экспериментальным значением 29 мНм~1. .
В [2] методом МД моделировали поверхность вода-нитробензол. Для обоих веществ были использованы гибкие эмпирические модели с парными потенциалами. Рассчитали профили плотности, ориентационные функции распределения, локальный электрический потенциал. Авторы [2] обнаружили, что хотя нитробензол является полярной жидкостью, структура поверхностного слоя на границе с водой аналогична структуре поверхностных слоев вода-органическая неполярная жидкость. Поверхностный слой характеризуется значительной шероховатостью с одновременно локально четкой границей. Дипольные моменты молекул воды предпочтительно расположены параллельно поверхности.
Работа [3] посвящена детальному изучению методом МД межфазной поверхности вода-закритическая двуокись углерода при давлениях 20 и 28 атм и температурах 318 и 328 К. Двуокись углерода представляет практический интерес как заменитель летучих органических жидкостей в процессах разделения и очистки веществ. Были рассчитаны профили плотности, корреляционные функции в поверхностном слое и объеме, ориентационные функции распределения, а также поверхностное натяжение. Последнее оказалось примерно на 50% выше экспериментального значения. На внутреннем профиле плотности наблюдали лики в поверхностном слое со стороны обеих жидкостей. В целом вывод о локальной структуре поверхностного слоя в [3] совпадает с результатом, полученным ранее для границы вода-масло [4], о том, что поверхностный слой между водой и нерастворимой жидкостью является достаточно узким с толщиной порядка нескольких ангстрем.
Гидрофильно-гидрофобная поверхность вода-гексан была рассмотрена в одной из недавних работ [5], где методом МД моделировали плоскую поверхность при температурах 283, 303, 323, 343 К. Рассчитанное поверхностное натяжение хорошо согласуется с экспериментальным значением.
Во всех указанных выше работах имели дело с плоскими поверхностями раздела двух жидких фаз. Однако часто одна из жидкостей находится в дисперсном состоянии, и тогда межфазный поверхностный слой обладает заметной кривизной. Поэтому настоящее исследование посвящено компьютерному моделированию кластеров воды,
© Г. В. Муджикова, Е. Н. Бродская, 2004
помещенных в неполярную жидкость. Свойства изолированных кластеров воды разного размера методом МД [6] подробно рассмотрены в работе [7] с использованием трех различных моделей воды БТ2 [8], Т1Р4Р [9] и БРС/Е [10]. Размеры кластеров варьировались в пределах от 64 до 1000 молекул воды. Наряду со структурными характеристиками оценивали работу образования кластеров воды из пара и определяли зависимость эффективного поверхностного натяжения от размера кластера. Было установлено, что поверхностные свойства воды очень чувствительны к выбору модели.
Задача данного исследования состоит в том, чтобы проследить, как влияет неполярная фаза на состояние кластера воды по сравнению с изолированным кластером.
В качестве неполярной фазы рассматривали жидкий н-гексан, как и в работе [5]. Все расчеты проводились методом МД в условиях 1ЧУТ ансамбля при 300 К для трех систем, в которых кластеры воды содержали 27, 64 и 125 молекул, в дальнейшем называемых системами I, II и III соответственно. Использовали алгоритм интегрирования уравнений движения с жесткими связями [6], временной шаг при этом составлял 1 фс. Вся система помещена в сферическую оболочку с отталкивательным потенциалом, радиус оболочки выбирался таким образом, чтобы средняя плотность гексана соответствовала плотности его жидкой фазы. Время численного эксперимента складывается из времени приведения системы к равновесию и времени усреднения. Первое составляло около 300 пс, второе - 2000 пс.
Для описания воды была применена трехточечная (БРС/Е) модель [10]. Она достаточно проста, что обеспечивает высокую эффективность вычислений, а объемную жидкость описывает так же хорошо, как и более сложные модели. Молекула представляет собой треугольную конструкцию из двух атомов водорода и атома кислорода. Длина валентной связи О-Н равна 0,1 нм, величина валентного угла НОН - 109,47°, связи жесткие, колебательные степени свободы не учитываются. Центры атомов в молекуле воды совпадают с положениями эффективных точечных зарядов. Помимо электростатического взаимодействия, представленного кулоновским потенциалом, между атомами кислорода действуют ван-дер-ваальсовские силы, описываемые эффективным парным потенциалом Леннард-Джонса (12-6). Молекулы гексана моделируются как простые сферические леннард-джонсовские частицы с параметрами, взятыми из работы [11, с. 390]. Взаимодействия частиц разного сорта рассчитывались при помощи правила Лоренца-Бертло.
Все используемые параметры для потенциалов приведены в табл." 1. Радиусы сферической оболочки составляли 2,9, 3,15 и 3,25 нм для систем I, II и III соответственно.
Таблица 1. Параметры моделей используемых потенциалов
Параметр qi, е Ei/кв, К CTj, нм
О -0,8476 78,172 0,3166
Н . 0,4238 - _
СбНн - 413 0,5909
Рассчитывались радиальные профили плотности воды pw(r) и гексана ph{r), профили кинетической и полной энергии е(г), ориентационные функции распределения молекул воды р(6) и профили локального нормального давления i^i(г). На основе последнего определялась работа образования кластера воды W и оценивалось эффективное поверхностное натяжение 7Эфф кластера воды в гексане. Также рассматривалась зависимость свойств системы от размера водного кластера, т.е. от числа молекул в нем. Ориентацию воды в микрокапле исследовали по данным для вектора дипольного момента молекулы и вектора, связывающего положения атомов водорода.
На рис. 1 показан радиальный профиль плотности вода-гексан для всех систем. Средние значения плотности в кластере воды оказываются выше по сравнению с экспериментальными в его объемной фазе (табл. 2). Это свидетельствует о том, что кластер
р
1,2 1,0 0,8
0,6 0,4 0,2 0,0
, г/см3
—1
—2
.у v
ч.' I
-1___I___ü-
0,0 0,4 0,8 1,2 1,6 2,0 2,4
Рис. 1. Профили плотности воды (1) и гек-сана (2) для систем I (о), II (б), III (в).
0,0 0,4 . 0,8 1Д 1,6 2,0 2,4
0,0 0,4 0,8 1,2 1,6 2,0 2,4
г, нм
воды испытывает большие давления со стороны фазы гексана. Хотя средняя плотность гексана близка значению в объеме, на профиле локальной плотности вблизи кластера воды проявляются значительные осцилляции. В области межфазной поверхности для всех систем прослеживается существенное взаимное пересечение профилей примерно на протяжении 0,7-0,8 нм, которое, однако, обусловлено не смешиваемостью компонентов, а значительными флуктуациями поверхности кластера за счет капиллярных волн. При этом ширина взаимного проникновения двух кривых характеризует степень отклонения поверхности кластера от сферической формы и соответственно амплитуду этих волн. На рис. 2 а, б представлены двухмерные проекции мгновенной конфигурации системы III, иллюстрирующие степень несферичности кластера.
Таблица 2. Данные по объемной, плотности для исследуемых подсистем
Плотность, г/см3 Система Эксперимент
I II III
Вода Гексан 1,108±0,017 0,573±0,005 1,008±0,031 0,559±0,006 1,070±0,008 б,535=Ь0,009 1,0 0,655
Помимо локальной плотности молекулярную структуру поверхностного слоя воды определяют ориентационные функции распределения р{в). Расположение молекул воды в поверхностном слое кластера можно охарактеризовать ориентацией вектора дипольного момента молекулы и вектора НН, проведенного через положения атомов водорода. Были рассчитаны функции р($) для вектора дипольного момента молекулы. У большинства молекул воды дипольный момент ориентирован нормально к радиус-
Рис. 2. Двухмерные проекции мгновенной конфигурации системы III.
вектору молекулы, при этом для вектора НН какой-либо преимущественной ориентации не наблюдается, что означает достаточно свободное вращение молекулы относительно своей оси. Для всех рассмотренных систем был получен одинаковый вид зависимости р{9). Таким образом, наиболее вероятное расположение молекул воды в поверхностном слое на границе с гидрофобной фазой остается таким же, как и на границе вода-пар: диполи воды ориентируются параллельно поверхности раздела.
При определении энергетических характеристик были рассмотрены профили локальной энергии исследуемых систем. Профили потенциальной энергии е(г) на молекулу демонстрирует рис. 3 для каждого компонента. Постоянство кинетической энергии (рис. 3, а, кривые 1 и 3) свидетельствует о термическом равновесии- в системе. Аналогичный результат был получен для всех систем. Потенциальная энергия у воды ниже, чем у гексана, что на качественном уровне хорошо согласуется с экспериментом, при этом минимальное значение для воды более чем в два раза ниже. Минимум е(г) для гексана находится в поверхностном слое, что может быть связано с более высокой плотностью его в области, примыкающей к поверхностному слою. Следует отметить, что в кластере воды и гексане отсутствуют области однородности по энергии, за исключением системы III.
На рис. 4 представлен профиль нормального давления для всех систем. При расчете этой характеристики использовали выражение для тензора давления Ирвинга-Кирквуда [7]. В отличие от нормального давления в изолированном кластере [7] в исследуемой системе отсутствует характерный для искривленного поверхностного слоя вода-пар переход в отрицательную область. Для всех систем наблюдается локальный максимум в межфазном поверхностном слое, что, по-видимому, можно объяснить сравнительно высоким внешним давлением в системах. Подтверждением этому является вид функций в области фазы гексана. Пренебрегая неоднородностями структуры гексана, на основе Р5ч оценивали среднее значение давления в нем. Для систем II и III оно составляет соответственно 1,12 • 107 и 0,8 ■ 107 Н/м2.
е-1019, Дж 0,0
-0,2
-0,4-, -0,6 -0,8 -1,0
» \ /' л \ /
: • \'
ь—
/
—1 --2 • • ■ 3 — 4
0,0 0,4 0,8 1,2 1,6 . 2,0 2,4
г, нм
Рис. 3. Профили парциальных энергий на молекулу для систем I (о), II (б) и III (в). 1,3- профили кинетической энергии молекул воды (1) и гексана (3) для системы I; 2, 4 ~ профили потенциальной энергии воды (2) и гексана
(4).
е-1019, Дж 0,0
-0,2
-0,4
-0,6
-0,8
-1,0
\ /
\ I _________
т-
0,0 -од
-0,4 -0,6 -0,8
0,0 0,4 0,8 1,2 1,6 2,0 2,4 в
-1,0
\ /
v ' \
— Г I
Ч-.__
0,0 0,4 0,8 1Д 1,6 2,0 2,4
г, нм
Ри-10"8, Нм""2 20
0,0 0,5 1,0 1,5 2,0 2,5 г, нм
Рис. 4• Профиль нормального давления (1) и его кинетическая составляющая ркТ (2) для систем I (а), II {б) и III (в).
Профиль давления позволяет рассчитать работу образования капли
Rß
W = 2тг J PN(r)r2dr - 2тгР^Л|/3, * (1)
о
где РР - давление в неполярной фазе; Rß - радиус некоторой сферы в области углеводородной фазы, для систем II и III он принят равным 2,63 и 2,73 нм, что соответствует работе 7,8 • Ю-20 и 15,6 • Ю~"20 Дж. Исходя из (1), по формуле
7эфф = ЗИ74тгД2 (2)
можно оценить эффективное поверхностное натяжение кластера для эквимолекулярной поверхности. Полученные результаты приведены в табл. 3 вместе с данными из [7]. Значения 7эфф для изолированного кластера и кластера в гексане оказались достаточно близкими. Как и в случае изолированного кластера, наблюдается тенденция роста эффективного поверхностного натяжения с увеличением размера кластера. Следует отметить, что полученные значения меньше, чем приведенное в работе [4] для плоской поверхности вода-гексан, которое составляет 48 мН/м.
Таблица 8. Оценка эффективного поверхностного натяжения и радиуса эквимолекулярной поверхности межфазной границы вода—н-гексан и вода-пар [6]
Параметр Вода-гексан Вода-пар
N (Н20) 64 .125 64 125
7эфф, мН/м 29,0 39,0 24,0 52,0
Re, HM 0,8 0,98 0,77 0,97
В результате численного эксперимента можно сделать следующие выводы. В отличие от изолированного кластера воды в изученных системах практически отсутствует взаимное проникновение молекул воды и гексана в соседнюю фазу. Кажущееся смешение в переходной зоне, наблюдаемое на радиальных профилях плотности, связано с флуктуациями формы кластера, характеризующие вклад капиллярных волн в толщину поверхностного слоя. По толщине переходной зоны вода-гексан на радиальных профилях плотности, полученных в предположении сферической симметрии системы, можно оценить среднюю величину флуктуаций, или амплитуду капиллярных волн.
Наиболее вероятное расположение молекул воды в поверхностном слое кластера в гексане такое же, как в изолированном кластере, т.е. диполи молекул ориентированы преимущественно вдоль поверхности. В поверхностном слое четко выделяется область повышенного давления, а также отсутствует свойственное системе вода-пар понижение до отрицательных значений. Это обусловлено достаточно высоким внешним давлением. Тенденция изменения эффективного поверхностного натяжения с увеличением размера кластера такая же, как и у кластера в вакууме [7], т. е. значение 7Эфф повышается с ростом кластера. При этом оно остается ниже, чем полученное в [5] поверхностное натяжение плоской поверхности.
Работа выполнена при финансовой поддержке программы Президента РФ «Ведущие научные школы» (грант № НШ-789.2003.3).
Summary
Mudzhikova G. V., Brodskaya E. N. Molecular dynamics simulation of small water clusters in nonpolar n-hexane.
Water clusters in liquid hexane were modeled by the molecular dynamics method at 300 K. The dependence of properties on cluster size was considered with clusters containing 27, 64 and 125 molecules. The radial profiles of the local density and the normal component of the pressure tensor were calculated. The effective surface tension was evaluated on the basis of the normal pressure profile. It was shown that there is a tendency to the increase of the surface tension along with the cluster size as in the case of the isolated water cluster.
Литература
1. Benjamin I. // J. Chem. Phys. 1992. Vol. 97. P. 1432-1445. 2. Michael D., Benjamin I. U J". Electroanal. Chem. 1998. Vol. 450. P. 335-345. 3. Da Rocha S. R. P., Johnston К. P, Westacott R. E., Rossky P. J. // J. Phys. Chem. B. 2001. Vol. 105. P. 12092-12104. 4. Senapati S., Berkowitz M. L. (I Phys. Rev. Lett. 2001. Vol. 87, N 17. P. 176101-176104. 5. Nicolas J. P., de Souza N. R. // J. Chem. Phys. 2004. Vol. 120. P. 2464-2469. 6. Allen M. P., Tildesley D. J. Computer simulation of liquids. Oxford, 1987. 7. Zakharov V. V., Brodskaya E. N., Laaksonen A. U J. Mol. Phys. 1998. Vol. 95. P. 203-209. 8. Stillinger F. H., Rahman A. // J. Chem. Phys. 1974. Vol. 64. P. 1545-1557. 9. Jorgensen W. L., Chandrasekhar J., Modura J. D. // J. Chem. Phys. 1983. Vol. 79. P. 926-935. 10. Berendsen H. J. C., Grigera J. R., Straatsma T. P. // J. Chem. Phys. 1987. Vol. 91. P. 6269-6271. 11. Справочник химика: В 10 т. 2-е изд. / Под ред. Б. П. Никольского. Л., 1962. Т. 1.
