Молекулярно-динамическое моделирование метастабильных фазовых состояний. Термодинамические свойства леннард-джонсовской системы Текст научной статьи по специальности «Физика»

ISSN 1992-6502 (P ri nt)_

2014. Т. 18, № 1 (62). С. 214-223

Ъыьмт QjrAQnQj

ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru

УДК 536.421

Молекулярно-динамическое моделирование метастабильных фазовых состояний. Термодинамические свойства леннард-джонсовской системы

1 2 3

С. П. Проценко , В. Г. Байдаков , З. Р. Козлова

1 [email protected], 2 baidakov@ itp.uran.ru , 3 [email protected] ФГБУН «Институт теплофизики Уральского отделения РАН» Поступила в редакцию 17.02.2014

Аннотация. Методом молекулярной динамики рассчитаны давление, внутренняя энергия и изохорная теплоемкость для 209 состояний леннард-джонсовской системы в интервале температур 0.35 < квТ/е < 2.0 и плотностей 0.001 < рст < 1.2. Полученный массив данных, наряду со стабильными состояниями, включает и однородные метастабильные (пересыщенный пар, перегретая и переохлажденная жидкость, перегретый кристалл). По данным моделирования построены термическое и калорическое уравнения состояния. Аппроксимированы спинодали пересыщенного пара и перегретой жидкости. В стабильной области полученные данные сопоставляются с результатами предшествующих работ. Для достижения требуемой точности рассчитываемых свойств усреднение микроскопических аналогов функций динамических переменных проводилось по 5-10 шагам интегрирования уравнений движения при моделировании конденсированных состояний и по (1-4)-10 шагам при исследовании газовой фазы. В однопроцессорном режиме моделирование одного состояния жидкости и кристалла составляло от 300 до 900 мин. Моделирование на системе из 64 процессоров при эффективности алгоритма и 0.5 сокращало время расчета примерно в 30 раз.

Ключевые слова: компьютерное моделирование; молекулярная динамика; термодинамические свойства; уравнение состояния.

ВВЕДЕНИЕ

Молекулярные системы, состоящие из большого числа частиц, могут находиться в нескольких фазовых состояниях. В определенной области температуры и давления фазы равновесно сосуществуют. При изменении состояния вдоль линии, пересекающей кривую равновесия фаз, вещество может остаться однородным, то есть сохранить свойства исходной фазы. Такие состояния не являются полностью устойчивыми. Они представляют метастабильные фазы вещества. Быстрые фазовые переходы в энергонапряженных процессах сопровождаются значительным отклонением системы от условия равновесия фаз. При этом одна из фаз оказывается в метастабильном состоянии. Чем интенсивнее процесс, тем больше отступление от по-

Работа выполнена в рамках программы Президиума РАН № 18 (проект 12-П-2-1049). Статья рекомендована к публикации программным комитетом международной научной конференции "Параллельные вычислительные технологии 2013".

ложения равновесия, обеспечивающее необходимую «движущую силу» превращения. В этом смысле метастабильные состояния не только возможны, но и неизбежны.

Информация о поведении вещества в мета-стабильных состояниях важна для разработки эффективных и безопасных технологических процессов в химии и энергетике, при получении, хранении и транспортировке криогенных жидкостей, создании новых конструкционных материалов, для описания кинетики фазовых и полиморфных переходов, кинетики разрушения материалов при высоких скоростях деформации. Актуальны исследования термодинамических и транспортных свойств в метастабильных состояниях, границ достижимого пересыщения метастабильных фаз, кинетики и молекулярных механизмов кристаллизации, плавления и кавитации, определение параметров фазового равновесия и поверхностного натяжения, механизмов зарождения и роста зародышей новых фаз.

Один из подходов к решению этих задач -натурный эксперимент. Экспериментальные исследования свойств и процессов в сильно ме-

тастабильных системах чрезвычайно затруднены малым временем их существования. Новые возможности здесь открывают методы компьютерного моделирования (методы Монте-Карло и молекулярной динамики). Они позволяют не только рассчитать свойства метастабильных фаз в экспериментально достижимых условиях, но и провести исследования фазовой метастабильно-сти при больших пересыщениях вблизи границы существенной неустойчивости - спинодали, а также исследовать кинетику фазовых превращений. Компьютерные модели позволяют реализовать степени метастабильности, значительно превышающие достигаемые в реальных экспериментах. Последнее обстоятельство обусловлено относительной малостью моделируемых систем. С другой стороны, с развитием многопроцессорных систем появилась возможность работать с моделями, содержащими десятки и сотни миллионов взаимодействующих частиц. Это позволяет приблизить результаты моделирования к свойствам макроскопических объектов, учесть корреляции на больших межчастичных расстояниях, присущие системам с сильной метастабильностью, описать свойства границ раздела равновесно сосуществующих фаз в плоском пределе. Решение всех указанных задач невозможно без реализации параллельных процессов на высокопроизводительных вычислительных системах. Новые перспективы перехода к исследованию систем больших масштабов открываются в связи с увеличением производительности вычислителей за счет оснащения GPU.

МОДЕЛЬ

Модели систем взаимодействующих частиц, основанные на описании взаимодействий потенциалом Леннард-Джонса (ЛД), традиционно широко используются для изучения простых флюидов. Потенциал ЛД учитывает наиболее существенные особенности поведения простого вещества не только в однородных состояниях (газ, жидкость, кристалл), но и пригоден для описания фазовых переходов и фазовых равновесий. Это свойство предопределяет его широкое использование для получения информации о свойствах простых веществ методами компьютерного моделирования с целью развития и совершенствования теорий жидкого состояния.

В данной работе моделирование проведено методом равновесной молекулярной динамики (МД) в NVE ансамбле, где N - число частиц, V - объем, E - внутренняя энергия системы. Все расчеты выполнены в кубической ячейке c

периодическими граничными условиями при N равном 2048 и 4000. Выбранные размеры систем обеспечивали возможность глубокого проникновения в область метастабильных состояний и их исследование за характерные времена компьютерных экспериментов. Частицы взаимодействовали посредством потенциала Леннард-Джонса, обрезанного при г = гс

Ф =

4 s

а | а

v r )

0,

r < r

r > r

(1)

где в и с - характеристические параметры энергии и расстояния.

Радиус обрезания потенциала выбирался равным г = 6.78с при плотностях р* = рс3 < 0.82 и был равен половине ребра ячейки для системы из 2048 частиц при больших плотностях, но не меньше чем 5.975с (р* =рс3 = 1.2). Выбор радиуса обрезания г* > 5.975 обеспечил: во-первых, пренебрежимо малый скачек сил при г = г* , во-вторых, максимально возможный в масштабах исследованных систем учет корреляций на больших межчастичных расстояниях, принципиально важный при моделировании ме-тастабильных и околокритических состояний.

Далее все величины представлены в приведенных единицах. Приведенные параметры обозначены знаком (*): расстояние г * = г / в, температура Т* = квТ / в, плотность р* =рс3, потенциальная энергия и * = и / в, внутренняя энергия е* = е/в, давление р* = рс3/в, изохорная теплоемкость с* = * / кв, где кв - постоянная Больцмана, т - масса частицы. Для интегрирования уравнений движения частиц использовался алгоритм Верле [1]. Шаг интегрирования по времени АI* = в /т / с = 0.0046 при температурах Т > 0.35. При Т = 0.35, АГ* = 0.0023 .

Расчеты проведены с применением пакета параллельного молекулярно-динамического моделирования LAMMPS [2] на вычислительных системах Института математики и механики УрО РАН и Межведомственного суперкомпьютерного центра РАН. В пакет LAMMPS встроен объект для расчета изохорной теплоемкости.

Для оценки качества программы использованы следующие показатели: ускорение = где N - время исполнения распараллеленной программы на N процессорах, t1 - время исполнения программы на одном процессоре; эффек-

<

тивность Е = определяющая среднюю долю времени выполнения параллельного алгоритма, в течение которого процессоры реально используются для решения задачи. Результаты расчетов ускорения и эффективности в интервале Ыр от 1 до 256 для системы из 4000 частиц в случаях моделирования фаз с плотностями, отличающимися на 2 порядка (жидкость, р = 0.75, газ, р = 0.0075), показаны на рис. 1.

40

30 -

С/)

20

10 -

1.0

- 0.8

- 0.6

ш

- 0.4

- 0.2

0.0

40 80

120 160 МР

200 240

Рис. 1. Зависимость ускорения (1, 2) и эффективности (3, 4) от количества процессоров на задачу и числовой плотности моделируемого состояния. Открытые значки отвечают плотности р = 0.0075, закрытые - р = 0.75

Резкое различие показателей ускорения и эффективности при моделировании систем разной плотности связано с различной загрузкой процессоров. При использованном в модели радиусе обрезания потенциала межчастичного взаимодействия в системе с плотностью

р = 0.75 среднее число соседей каждой частицы

составляло 605 атомов, а при р = 0.0075 уменьшалось до 6. В соответствии с результатами вычислительного эксперимента, представленными на рис. 1, для моделирований конденсированных состояний использовалось 64 процессора, что обеспечивало ускорение, равное 30 при эффективности и 0.5 и продолжительности расчета до 40 мин. Моделирование газовой фазы проводилось не более чем на 32 процессорах с ускорением до 5 и эффективностью и 0.18.

МЕТОДИКА МОДЕЛИРОВАНИЯ МЕТАСТАБИЛЬНЫХ СОСТОЯНИЙ

Расчеты проводились по изотермам в интервале температур от 0.1 до 2.0 и плотностей от 0.001 до 1.2. В газе начальные состояния отвечали приведенной плотности 0.001, в ГЦК кристалле плотности 1.2. В жидкости на изотермах от 0.7 до 2.0 расчеты начинались из стабильной области состояний. При заданной температуре последовательный переход к следующим состояниям осуществлялся линейным масштабированием ребер ячейки (сжатием, растяжением) и координат частиц конечной конфигурации предыдущего состояния. По мере роста степени метаста-бильности шаг по плотности уменьшался. Процедура последовательного увеличения (уменьшения) плотности при постоянной температуре обеспечивала глубокие заходы за линии равновесия фаз без потери однородности системы. Резкое изменение состояния приводило к возникновению локальной неоднородности и преждевременному фазовому переходу. Время уравновешивания варьировалось от 100 до 250 тыс. шагов интегрирования уравнений движения и зависело от термодинамического состояния системы.

При температурах Т < 0.7 все состояния жидкой фазы метастабильны. Большинство из них относится к области отрицательных давлений. В качестве начальной конфигурации частиц для этих температур использовалась конфигурация, полученная при Т = 0.7 и р = 0.875. Переход к требуемой температуре осуществлялся последовательным изохорическим охлаждением. Масштабированием размеров ячейки и межчастичных расстояний создавались состояния с меньшей и большей плотностью.

Расчеты проводились до плотностей, при которых однородные состояния разрушались в результате спонтанного образования и последующего роста зародыша новой фазы [3, 4]. Распад метастабильной фазы фиксировался по скачку давления (см. рис. 2). В этот момент расчеты всех термодинамических параметров прекращались. Если до распада метастабильного состояния набранной статистики было недостаточно, то генерировалась новая МД траектория и расчеты повторялись. Расчет проводился до тех пор, пока время жизни метастабильного состояния не оказывалось достаточным для оценки термодинамических параметров с требуемой точностью.

Метастабильная фаза допускает наличие в ней кластеров новой фазы, размеры которых не превышают критического. Несмотря на присутствие кластеров докритических размеров, уро-

0

0

вень флуктуаций температуры, давления, потенциальной энергии и других термодинамических свойств в метастабильном состоянии незначительно отличается от их значений в области стабильных состояний. Появление зародыша закритического размера сразу сопровождается его необратимым ростом. Такие события легко фиксируются по скачкообразному изменению термодинамических свойств системы. Наиболее чувствительным свойством к образованию за-критических зародышей является давление. Характерное время разрушения метастабильного состояния с образованием двухфазной микрогетерогенной системы составляет всего несколько тысяч шагов интегрирования уравнений движения частиц.

-1

2 Ар =0.02

—' 3 Ар =0.016

1 JU. ll! . 1 ..Ь ^Л iu.il. ^ ..1 и 1

■ г. ' г ' гг' щмш Т|1 1Т 1 'п» И 1

■............"'" 4 1 Ар =0.02 1 1

0

500

1000

1500

Рис. 2. Зависимость давления от времени при Т = 0.85: 1 - р * = 1.015, переход жидкость-кристалл; 2 - р * = 1.01, переохлажденная жидкость; 3 - р = 0.85, стабильная жидкость; 4 - р * = 0.69 , перегретая жидкость; 5 - р * = 0.67, переход жидкость-пар. Прямыми горизонтальными

линиями показаны средние по времени

. *

значения давления, Ар - стандартное отклонение давления в однородных состояниях

На рис. 2 представлена зависимость давления от времени для Т = 0.85 при пяти значениях плотности. Рисунок демонстрирует постоянство (с точностью до флуктуаций) давления как в состоянии стабильной жидкости (кривая 3), так и в метастабильных состояниях жидкой фазы вблизи границ достижимого перегрева (кривая

4) и переохлаждения (кривая 2). Стандартные отклонения давления от средних значений в состояниях стабильной, перегретой и переохлажденной жидкости практически совпадают. Появлению и последующему росту как критического кристаллика (кривая 1), так и критического пузырька пара (кривая 5) отвечают скачки давления.

Поскольку амплитуда таких скачков в несколько раз превышает уровень флуктуаций в однофазных состояниях, моменты разрушения метастабильных состояний однозначно идентифицируются. Подобный контроль «однофазно-сти» входил в стандартную процедуру наших компьютерных экспериментов и проводился во всех исследованных состояниях. Это позволило при расчете термодинамических и транспортных свойств исключать двухфазные микрогетерогенные состояния. Использованный подход не противоречит результатам более детального анализа структуры систем до момента распада метастабильного состояния, основанного на выявлении кластеров новой фазы, определении их размеров и построении функции распределения кластеров по размерам. Этот метод был использован нами при исследовании кинетики заро-дышеобразования в метастабильных системах

[3, 4].

РЕЗУЛЬТАТЫ РАСЧЕТОВ ТЕРМОДИНАМИЧЕСКИХ СВОЙСТВ

В молекулярно-динамических экспериментах рассчитаны кинетическая, потенциальная и внутренняя энергии, температура, давление и изохорная теплоемкость в стабильных (газ, жидкость, кристалл) и метастабильных (пересыщенный газ, перегретая и перохлажденная жидкость, перегретый кристалл) состояниях леннард-джонсовской системы.

Кинетическая энергия на частицу равна

ек =

1 Ы — £

2 N £

*2

(2)

где V* - скорость 7-й частицы.

Потенциальная энергия на частицу есть

и =

■>=2 г* <г'*

8 лр

3г*

*3 •

N 1=1

Внутренняя энергия на частицы равна

(3)

е = е + и .

Температура рассчитывается согласно

7

6

1

0

Т = /3.

Давление определяется как

г

* р р =т

1 N-1 N 2е* --1 £ £ * N £ £

аф( г*)

1=2 . " ^

1 >', Гц йГ* у

*2

16 лр

Изохорная теплоемкость с* рассчитывалась через среднеквадратичные флуктуации кинетической энергии [5]

С /

*2 * 2 1 -1N <>-<^ >

Л-1

< е* >

(7)

Вклады в давление и потенциальную энергию частиц, находящихся на расстояниях

г > г* , учитывались поправками, которые вводились в приближении отсутствия корреляции во взаимодействии частиц на таких расстояниях. Вклад отталкивательного члена потенциала взаимодействия в поправках не учитывался. Для изохорной теплоемкости поправка на дальнодействие не вводилась. При плотностях, близких к критической плотности, вклад вторых слагаемых уравнений (3) и (6) в давление составляет « 4%, в потенциальную энергию 2.8%, а при максимальной плотности р* = 1.2: 0.3 и 1.2%

соответственно.

Время выхода системы на равновесие и время, по которому определялись средние значения термодинамических свойств, зависело от термодинамического состояния системы (стабильное, метастабильное, околокритическое), плотности, температуры, рассчитываемого свойства. В стабильных состояниях и в состояниях, удаленных от критической точки, для определения средних значений рассчитываемых величин с погрешностью «1% достаточно усреднения по десяткам тысяч состояний. В области сильной метастабильности такая точность достигается только увеличением интервала усреднения. Особенно критичной к интервалу усреднения была изохорная теплоемкость с*. В состояниях, близких к границам достижимых пересыщений, регистрировались выбросы частичных средних значений, полученных по ин-

о *

тервалам времени 5 • 10 А £. Результаты для с* с

(5) погрешностью « 1% во всем исследованном интервале плотностей и температур получены усреднением по (0.1 - 1)-106 микросостояний. В

*

ряде случаев получить надежные значения с при температурах ниже 0.5 не удалось вследствие распада метастабильной жидкости за харак-

(6) терные времена компьютерных экспериментов. С учетом этих обстоятельств, длина интервала уравновешивания выбиралась (100 - 250)-103 временных шагов, средние значения определялись по интервалам 5-105 Д£* для жидкости и 1-106 Д* для газа. При исследовании разреженных систем особое внимание необходимо уделять накоплению достаточной статистики столкновений. Поэтому все термодинамические свойства получены усреднением не менее чем по 106 состояний. Дополнительные пробные расчеты длиной 2-106 и 4-106 шагов показали совпадение средних значений любых термодинамических свойств. Таким образом, стандартной траектории длиной в 106 шагов достаточно для оценки свойств в газовой фазе с погрешностью не хуже чем в плотной жидкости.

Поскольку компьютерные эксперименты проводились в условиях постоянства N V, Е, температура не являлась фиксированным параметром. Поэтому при каждом значении плотности расчеты проводились для двух температур,

близких к заданной. Значения р, е, Су , соответствующие требуемой температуре, определялись линейной интерполяцией на температурных интервалах, не превышающих 0.03.

Погрешности средних значений р, и, с* рассчитывались посредством разбиения полных наборов мгновенных значений указанных величин на блоки длиной 5-103 значений. Частичные средние по блокам формировали массивы независимых измерений, необходимых для вычисления среднеквадратичных отклонений.

Результаты расчета давления, внутренней энергии и изохорной теплоемкости по изотермам представлены на рис. 3 - 8. Среднеквадратичные погрешности рассчитанных свойств укладываются в размеры значков на соответствующих графиках. Полученные данные по давлению и внутренней энергии в стабильных и метастабильных состояниях леннард-джонсовской системы использованы для построения локальных термических и калорических уравнений состояния

р =££«^т *

1=0 i=0

*

с

=££*, р* т4. (9)

1=0 i=0

Уравнение (8) позволяет оценить положения ветвей спинодали - границы существенной неустойчивости метастабильной фазы, определяемой условием

Г ^ * \ op

v°P* ^.

= о

(10)

0.20

0.15 -

0.10 -

0.05 -

0.00

0.0

0.1

0.2

0.3

P

Рис. 3. Изотермы давления (газ): 1 - Т = 0.55, 2 - 0.7, 3 - 0.85, 4 - 1.0, 5 - 1.15, 6 - 1.2, 7 - 1.25, 8 - 1.3, 9 - 1.5, 10 - 2.0. Сплошные линии - расчет по уравнению состояния (8) для газовой фазы. Штриховая линия - линия фазового равновесия газ-жидкость [6], С - критическая точка, СВ - спинодаль газа

Коэффициенты уравнений (8), а также максимальные значения показателей степеней п и т определялись методом регрессионного анализа. Для газовой фазы значения коэффициентов

равны

af2 = -13.9024,

af = -91.6712,

a0g4 = 525.696, a0g5 = -255.649, ag = 0.994318, ag = 16.0688, ag = 151.379, ag = -854.9177,

12 ,0

"13 i0 -

ag = 291.499, ag = 0.005776, af2 = -9.29945,

ag3 = 13.9024.

Для

al02 =-31.7465 43 = 75.0892,

жидкои

l

фазы:

a'os = 14.9614, a[4 =-61.2481,

alm = 12.727, a[2 =-18.1082, a[5 = 19.1732,

a2l = 1.80406, al22 =-5.71191, a\2 = 0.64104.

Для кристаллической фазы: a2° = -69.8834 a2 = 286.813, a^ = -412.723 ~cr a2cr0 = -31.8944, =-17.2507, a,°° = 71.9482, a'0 = -170.891

a;0 = 9.4696

a2° = 191.151, a°°= 47.2784,

a

22

a,°° = 135.637, a2° = -36.0205.

4 -

1 -

40

e - i

□ - 2 - 3

ф - 4

О - 5

< - 6

V - 7

О - 8

В - 9

О - 10

О - 11

■Й- - 12

О - 13

Ф - 14

Д - 1

- 35

30

25

20

a

15

10

- 5

P

Рис. 4. Изотермы давления (жидкость): 1 - T* = 0.35, 2 - 0.4, 3 - 0.45, 4 - 0.5, 5 - 0.55, 6 - 0.7, 7 - 0.85, 8 - 1.0, 9 - 1.15, 10 - 1.2, 11 - 1.25, 12 - 1.3, 13 - 1.35,

14 - 1.5, 15 - 2.0. Сплошные линии - расчет по уравнению состояния (8) для жидкой фазы, штриховые линии: линии фазового равновесия жидкость-газ [6] и жидкость-кристалл [7], штрих-пунктирная линия: спинодаль перегретой жидкости

Для сопоставления результатов наших расчетов с данными других авторов мы использовали уравнение состояния Kolafa and Nesbeda [8]. Данное уравнение построено относительно свободной энергии Гельмгольца и позволяет описывать как термические, так и калорические свойства леннард-джонсовского флюида. Область действия уравнения по температуре охватывает интервал от T = 0.68 до 10.

Уравнение состояния Kolafa and Nesbeda [8] качественно правильно воспроизводит поведение изотерм давления газа и жидкости как в стабильной, так и в метастабильной областях (см. рис. 3, 4). Хорошее количественное согласие по давлению имеет место для газовых ветвей изотерм, исключая участки, прилегающие к спинодали. Параметры спинодали газа, рассчи-

5

3

2

0

0

0.4

0.6

0.8

0.9

.0

.1

.2

*

танные из уравнения состояния Kolafa and Nesbeda [8] и локальных уравнений состояния (8), существенно расходятся. Так, по плотности на спинодали расхождения составляют ±(3 — 8)%. Аналогично со стороны жидкой фазы. Здесь вблизи спинодали расхождения по плотности составляют ±(1 — 3)%. Погрешность описания наших данных уравнением состояния [8] также растет по мере приближения к нижней температурной границе уравнения (T = 0.68) и с заходом в область переохлажденных состояний

жидкой фазы. При T* = 0.4 расхождение по дав* * *

лению Ap = pMD — pcalc с уравнением состояния [8] достигает -0.3, т. е. примерно (15 — 20)%. В области переохлажденных состояний жидкости эта величина достигает -0.6 при T * = 2.

20 -

15 -

10 -

5 -

❖ - 1

В - 2

в - 3 > - 4

О - 5

Д - 6

V - 7

О - 8

О - 9

if -10

□ -11

0.8

0.9

1.0

*

Р

1.1

1.2

Рис. 5. Изотермы давления (кристалл):

I - Т = 0.1, 2 - 0.2, 3 - 0.3, 4 - 0.4, 5 - 0.55, 6 - 0.7, 7 - 0.85, 8 - 1.0, 9 - 1.15, 10 - 1.5,

II - 2.0. Штриховая линия - линия фазового

равновесия жидкость-кристалл [7], ЛБ - спинодаль кристалла

Для внутренней энергии наибольшее расхождение с уравнением состояния работы [8] также наблюдается на метастабильных участках изотерм и при Т<0.68. Причем, в отличие от термического уравнения состояния, более существенные расхождения имеют место не столько в абсолютных значениях внутренней энергии, сколько в величине производной (де / др)г.

Причиной этих расхождений является то, что при построении уравнения состояния [8] в набор данных не были включены результаты рас-

четов в метастабильных состояниях леннард-джонсовской системы.

-2 -

ш

-4 -

-6 -

0.0

0.2

0.4

0.6

0.8

1.0

1.2

Р

Рис. 6. Изотермы внутренней энергии газа и жидкости: 1 - Т = 0.4, 2 - 0.55, 3 - 0.7, 4 - 0.85, 5 - 1.0, 6 - 1.15, 7 - 1.2, 8 - 1.25, 9 - 1.3, 10 - 1.35, 11 - 1.5, 12 - 2. Сплошные линии - результаты расчета по уравнению состояния [8], штриховые линии - линии фазового равновесия жидкость-газ [6]

-2 -

-4 -

ф

-6

О - 1

" в - 2

е - 3

- О ► - 4

о • - 5

. А ▲ - 6

V Т - 7

О ♦ - 8

" < < - 9

•к ★ - 10

- □ 1 ■ - 11 I

0.6 0.7

_L

_L

_L

_L

0.8 0.9 1.0 1.1 Р*

1.2

Рис. 7. Изотермы внутренней жидкости (темные точки) и кристалла (светлые точки): 1 - Т = 0. 1, 2 - 0.2, 3 - 0.3, 4 - 0.4, 5 - 0.55, 6 - 0.7, 7 - 0.85, 8 - 1.0, 9 - 1.15, 10 - 1.5, 11 - 2.0. Сплошные линии - результаты расчета по уравнениям состояния [8], штриховые линии - линии фазового равновесия жидкость-кристалл [7]

2

0

0

0

Результаты расчетов изохорной теплоемкости представлены на рис. 8. В отличие от механической, термическая устойчивость флюидной фазы в районе фазового перехода ведет себя более сложным образом (см. рис. 8). По мере захода в область метастабильных состояний механическая устойчивость жидкости и пара понижается (производная (ф/Зр) уменьшается).

3.2 3.0 2.8 2.6 2.4 2.2 2.0

1.8

1.6

* ^ О

0.0

0.3

0.6

0.9

1.2

Р

*

Рис. 8. Изотермы изохорной теплоемкости с* :

I - Т* = 0.4, 2 - 0.55, 3 - 0.7, 4 - 0.85, 5 - 1.0, 6 - 1.15, 7 - 1.2, 8 - 1.25, 9 - 1.3, 10 - 1.35,

II - 1.5, 12 - 2.0. Обозначения символов см. рис. 6. Штриховые линии - бинодаль [6],

штрих-пунктирные - спинодаль

Из рис. 8 следует, что если изохорная теплоемкость газовой фазы монотонно возрастает с ростом пересыщения (T = const), то теплоемкость жидкости по мере приближения к линии насыщения и заходом в метастабильную область проходит через минимум и далее возрастает при движении к спинодали.

ЗАКЛЮЧЕНИЕ

Леннард-джонсовский флюид является «пробным камнем» для теоретических моделей простых жидкостей. Это предопределяет актуальность изучения его термодинамических, кинетических и структурных характеристик в компьютерных экспериментах методами Монте-Карло и молекулярной динамики. Особый интерес в поведении молекулярной системы представляют точки фазовых переходов. Положение данных точек на термодинамической поверхно-

сти состояний определяется силами межмолекулярного взаимодействия. Универсальность микроскопического поведения различных веществ при сжатии от идеально-газового состояния с превращениями сначала в жидкость, а затем в кристалл, отражает универсальность взаимодействия частиц при изменении расстояния между ними.

Нами методом молекулярной динамики рассчитаны давление, внутренняя энергия газа, жидкости и кристалла и изохорная теплоемкость в области фазового перехода жидкость-газ. Как любой переход первого рода, фазовый переход жидкость-газ сопровождается фазовой метастабильностью, а именно - перегревом жидкости и переохлаждением газовой фазы. Так как вероятность появления в метастабильной системе за заданный промежуток времени зародыша новой фазы прямо пропорциональна ее объему, то малые размеры молекулярно-динамической системы позволяют осуществить глубокий заход в метастабильную область и рассчитать здесь термодинамические и другие свойства. В серии компьютерных экспериментов получены данные о термодинамическом потенциале, его первой и второй производной для 209 состояний, из которых 101 относятся к области метастабильных состояний. Самая низкая температура, при которой удалось удержать ме-тастабильную жидкую фазу без нарушения ее однородности, составила Т = 0.35 . При данной температуре жидкость находится в растянутом состоянии при отрицательном давлении

р =-2.3.

Характер плотностной зависимости давления на изотермах подобен предсказываемому уравнением состояния ван-дер-Ваальса. Однако вторая производная термодинамического потенциала (изохорная теплоемкость) свидетельствует о принципиальном отличии в термодинамическом поведении ван-дер-ваальсовского и леннард-джонсовского флюидов. Изохоры давления ван-дер-ваальсовского флюида прямолинейны, а изохорная теплоемкость не зависит от плотности, в то время как изохоры леннард-джонсовского флюида имеют кривизну, которая меняет свой знак по мере приближения к линии фазового перехода и заходом в метастабильную область. О последнем убедительно свидетельствуют изотермы изохорной теплоемкости, которые, как видно из рис. 8, содержат минимумы, максимумы и точки перегиба. Такой характер поведения изохорной теплоемкости в районе фазового перехода жидкость-газ качественно

*

и количественно согласуется с результатами прямых измерений в аргоне.

Имеющиеся в литературе уравнения состояния леннард-джонсовского флюида, построенные по результатам компьютерных экспериментов, качественно правильно передают зависимость давления от плотности и температуры не только в стабильной, но и в метаста-бильной области. Однако количественные расхождения с результатами наших расчетов увеличиваются по мере роста пересыщения мета-стабильной фазы и приближения к спинодали, а также при понижении температуры, так как во всех предшествующих работах нижний температурный предел был ограничен температурой тройной точки. Более существенные рассогласования эмпирических уравнений состояния и наших данных имеют место для изохорной теплоемкости. Здесь расхождения достигают 50 % и более. Тем не менее, среди эмпирических уравнений состояния леннард-джонсовского флюида можно выделить уравнения состояния [8, 9], которые наиболее лучшим образом согласуются с нашими данными по первым и вторым производным термодинамического потенциала.

Полученные в ходе компьютерного эксперимента термодинамические свойства леннард-джонсовского флюида использованы для составления локальных уравнений состояния. Такие уравнения, описывающие стабильные и ме-тастабильные области, необходимы для определения положения спинодали, описания кинетики нуклеации, проверки микроскопических теорий фазовой метастабильности.

СПИСОК ЛИТЕРАТУРЫ

1. Verlet L. Computer «Experiments» on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules // Phys. Rev. 1967. Vol. 159, No. 1. P. 98-103. [ L. Verlet, "Computer «Experiments» on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules," Phys. Rev., vol. 159, no. 1, pp. 98-103, 1967. ]

2. LAMMPS Molecular Dynamics Simulator [Электронный ресурс]. URL: http://lammps.sandia.gov (дата обращения 05.10.2012). [ (2012, Oct. 05) LAMMPS Molecular Dynamics Simulator [Online]. Available: http://lammps.sandia.gov ]

3. Protsenko S. P., Baidakov V. G., Teterin A. S., Zhdanov E. R. Computer simulation of nucleation in a gas-saturated liquid // J. Chem. Phys. 2007. Vol. 126, No. 9. P. 094502(1-14). [ S. P. Protsenko, V. G. Baidakov, A. S. Teterin, E. R. Zhdanov, "Computer simulation of nucleation in a gas-saturated liquid," J. Chem. Phys., vol. 126, no. 9, pp. 094502(1-14), 2007. ]

4. Baidakov V. G., Tipeev A. O., Bobrov K. S., Ionov G. V. Crystal nucleation rate isotherms in Lennard-Jones liquids // J. Chem. Phys. 2010. Vol. 132, No. 23. P. 234505 (1-9). [ V. G. Baidakov, A. O. Tipeev, K. S. Bobrov, G. V. Ionov, "Crystal nucleation rate isotherms in Lennard-Jones liquids," J. Chem. Phys., vol. 132, no. 23, pp. 234505 (1-9), 2010. ]

5. Lebowitz J. L., Percus J. K., Verlet L. Ensemble Dependence of Fluctuations with Application to Machine Computations // Phys. Rev. 1967. Vol. 153, No. 1. P. 250-254. [ J. L. Lebowitz, J. K. Percus, L. Verlet, "Ensemble Dependence of Fluctuations with Application to Machine Computations," Phys. Rev., vol. 153, no. 1, pp. 250-254, 1967. ]

6. Baidakov V. G., Protsenko S. P., Kozlova Z. R., Chernykh G. G. Metastable extension of the liquid-vapor phase equilibrium curve and surface tension // J. Chem. Phys. 2007. Vol. 126, No. 21. P. 214505(1-9). [ V. G. Baidakov, S. P. Protsenko, Z. R. Kozlova, G. G. Chernykh, "Metastable extension of the liquid-vapor phase equilibrium curve and surface tension," J. Chem. Phys., vol. 126, no. 21, pp. 214505(1-9), 2007. ]

7. Baidakov V. G., Protsenko S. P. Singular Point of a System of Lennard-Jones Particles at Negative Pressures // Phys. Rev. Lett. 2005. Vol. 95, No. 1. P. 015701(1-4). [ V. G. Baidakov, S. P. Protsenko, "Singular Point of a System of Lennard-Jones Particles at Negative Pressures," Phys. Rev. Lett., vol. 95, no. 1, pp. 015701(1-4), 2005. ]

8. Kolafa J., Nesbeda I. The Lennard-Jones fluid: an accurate analytic and theoretically-based equation of state // Fluid Phase Equilibria. 1994. Vol. 100. P. 1-34. [ J. Kolafa, I. Nesbeda, "The Lennard-Jones fluid: an accurate analytic and theoretically-based equation of state," Fluid Phase Equilibria, vol. 100, pp. 1-34, 1994. ]

9. Mecke M., Müller A., Winkelmann J., Vrabec J., Fischer J., Span R., Wagner W. An Accurate Van der Waals-Type Equation of State for the Lennard-Jones Fluid // Int. J. Thermophysics. 1996. Vol. 17, No. 2. P. 391-404. [ M. Mecke, A. Müller, J. Winkelmann, J. Vrabec, J. Fischer, R. Span, W. Wagner, "An Accurate Van der Waals-Type Equation of State for the Lennard-Jones Fluid," Int. J. Thermophysics, vol. 17, no. 2, pp. 391-404, 1996. ]

ОБ АВТОРАХ

ПРОЦЕНКО Сергей Павлович, ст. науч. сотр. лаб. криогени-ки и энергетики. Дипл. инж.-физик (УПИ, 1972). Канд. физ.-мат. наук (УПИ, 1978). Иссл. в обл. комп. моделирования молекулярных систем.

БАЙДАКОВ Владимир Георгиевич, зав. лаб. криогеники и энергетики. Дипл. инж.-физик (УПИ, 1970). Д-р физ.-мат. наук (ФТИНТ, 1987). Иссл. в обл. метастабильных состояний и комп. моделирования молекулярных систем.

КОЗЛОВА Залия Рафиковна, мл. науч. сотр. лаб. криогеники и энергетики. Дипл. учитель физики и информатики (Курганск. гос. ун-т, 1998). Иссл. в обл. комп. моделирования молекулярных систем .

METADATA

Title: Molecular-dynamics simulation of metastable phase states. Thermodynamic properties of Lennard-Jones system.

Authors: S. P. Protsenko1, V. G. Baidakov2, Z. R. Kozlova3 Affiliation: Institute of Thermophysics, Ural Branch of the

Russian Academy of Sciences, Russia. Email: 1 [email protected]. Language: Russian.

Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 18, no. 1 (62), pp. 214-223, 2014. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print).

Abstract: The method of molecular dynamics has been used to calculate the pressure, internal energy and isochoric heat capacity of the Lennard-Jones fluid for 209 states in the range of temperatures 0.35 < kBT/g < 2.0 and densities 0.001 < pa3 < 1.2. The array of data obtained, along with stable states, includes homogeneous metastable states (supersaturated vapor, superheated and supercooled liquid). Spinodals of supersaturated vapor and superheated liquid have been approximated. In a stable region the data obtained are compared with the results of previous papers.

Key words: Computer simulation; molecular dynamics; thermodynamic properties; equation of state.

About authors:

PROTSENKO, Sergey Pavlovich, Senior scientist, cryogenics and energetics laboratory. Dipl. engineer-physicist (Ural Politechnical inst. 1972). Cand. of Phys.-Math. Sci. (UPI, 1978).

BAIDAKOV, Vladimir Georgievich, Head of Laboratory cryogenics and energetics. Dipl. engineer-physicist (Ural Politechnical inst. 1970). Dr. of Phys.-Math. Sci. (FTINT, 1987).

KOZLOVA, Zaliya Rafikovna, Junior Researcher, cryogenics and energetics laboratory. Dipl. teacher of physics and informatics (Kurgan State Univ., 1998).