1. На главную
  2. Электронные версии
  3. Теплофизика высоких температур
  4. T. 57, № 1, 2019

Теплофизика высоких температур, 2019, T. 57, № 1, стр. 32-36

Расчет сжимаемости аргона при различных скоростях охлаждения

Е. И. Герман 1,  2, *, Ш. Б. Цыдыпов 1, **, Б. Б. Дамдинов 1

1 Бурятский государственный университет,
г. Улан-Удэ,, Россия

2 Геологический институт СО РАН,
г. Улан-Удэ,, Россия

* E-mail: net-admin@list.ru
** E-mail: shulun@bsu.ru

Поступила в редакцию 31.10.2017
После доработки 10.10.2018
Принята к публикации 13.07.2018

Полный текст (PDF)

Аннотация

Предложены приемы моделирования термодинамических процессов в молекулярных системах. Приведены результаты расчета изотермической сжимаемости аргона методом молекулярной динамики при изобарном охлаждении от 200 до 15 K при давлении 40 атм, адиабатической сжимаемости в диапазоне от 200 до 160 К при 40 атм. На температурных зависимостях изотермической сжимаемости при быстром (1012 К/с) и относительно медленном (109 К/с) охлаждении видны скачки, трактуемые как фазовые переходы.

ВВЕДЕНИЕ

Теплофизические и другие свойства стеклообразных твердых тел во многом определяются структурой жидкости, которая фиксируется в области стеклования. Вместе с тем в настоящее время нет общепризнанной теории перехода жидкость–стекло, а также нет полной ясности в природе жидкого и стеклообразного состояний вещества. Построение строгой статистической теории встречает большие трудности. Успехи, достигнутые в последнее время, касаются лишь равновесных свойств простых жидкостей со сферически симметричным парным потенциалом взаимодействия. Поэтому на данном этапе остается актуальным развитие моделей структуры жидкостей и стекол, основанных на численных и натурных экспериментах [1–3].

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

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

МЕТОДИКА МОДЕЛИРОВАНИЯ

Методика МД основана на численном интегрировании дифференциальных уравнений движения для системы многих частиц. Сила, действующая на i-ю частицу, связана с потенциалом взаимодействия Φ(r):

$m\frac{{{{d}^{2}}{\mathbf{r}}(t)}}{{d{{t}^{2}}}} = - \sum\limits_j^{N - 1} {\frac{{d{\text{Ф }}({{{\mathbf{r}}}_{{ij}}})}}{{d{\mathbf{r}}}}} ,$
где m – масса частицы, r – радиус-вектор частицы. Для системы инертных частиц хорошо подходит модель взаимодействия, предложенная Леннардом-Джонсом:
${\text{Ф }}({{{\mathbf{r}}}_{{ij}}}) = 4\varepsilon \left( {{{{\left( {\frac{\sigma }{{{{r}_{{ij}}}}}} \right)}}^{{12}}} - {{{\left( {\frac{\sigma }{{{{r}_{{ij}}}}}} \right)}}^{6}}} \right),$
где ε – глубина потенциальной ямы, σ – эффективный размер частицы.

Для численного решения уравнений движения хорошо себя зарекомендовал алгоритм Верле в скоростной форме, согласно которому положение частиц и их скорости на каждой итерации вычисляются в следующей последовательности:

${{{\mathbf{r}}}_{i}}(t + \Delta t) = {{{\mathbf{r}}}_{i}}(t) + {{v}_{i}}(t)\Delta t + {{{\mathbf{a}}}_{i}}(t)\frac{{\Delta {{t}^{2}}}}{2},$
${{{\mathbf{v}}}_{i}}\left( {t + \frac{{\Delta t}}{2}} \right) = {{{\mathbf{v}}}_{i}}(t) + \frac{1}{2}{{{\mathbf{a}}}_{i}}(t)\Delta t,$
${{{\mathbf{a}}}_{i}}(t + \Delta t) = - \frac{1}{m}\sum\limits_j^{N - 1} {\frac{{d{\text{Ф }}({{{\mathbf{r}}}_{{ij}}})}}{{d{\mathbf{r}}}}} ,$
${{{\mathbf{v}}}_{i}}(t + \Delta t) = {{{\mathbf{v}}}_{i}}\left( {t + \frac{{\Delta t}}{2}} \right) + {{{\mathbf{a}}}_{i}}(t + \Delta t)\frac{{\Delta t}}{2},$
где ri, vi, ai – радиус-вектор положения, скорость и ускорение i-й частицы соответственно.

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

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

Другой метод неравновесной динамики, основанный на введении возмущений на этапе инициализации системы, описан в [7, 8]. Ключевым моментом данной методики являются изучение релаксационных свойств системы и интерпретация полученных результатов для описания свойств сплошной среды.

Для определения закономерностей, происходящих в точках фазовых переходов в [9, 10], произведены компьютерные эксперименты с инициализацией молекулярной системы в области межфазового перехода. Подобным же образом исследовалась структура “замороженных” систем (стекол) в работах [11–13]. Стоит отметить, что структура систем, инициированных таким образом в условиях твердой фазы несопоставима со структурой реальных физических систем, так как в них образование структуры кристалла или аморфного твердого тела происходит в процессе остывания жидкости в зависимости от скорости охлаждения.

Компьютерные эксперименты с охлаждением жидкости и последующей кристаллизацией (в случае медленного охлаждения) или стеклованием (быстрое охлаждение) и расчетами взаимосвязи их упругих характеристик приведены в [14–19]. В указанных работах моделирование охлаждения системы частиц в широком температурном диапазоне производится при постоянном объеме системы. Это связано с тем, что методика моделирования изохорного охлаждения предполагает довольно простое изменение классического алгоритма Верле вводом процедуры пропорционального изменения всех скоростей частиц системы через заданное количество итераций. Такой искусственный способ контроля кинетической энергии системы критикуется в [5], где предлагается моделировать взаимодействие системы с термостатом через поверхность (использование стохастических граничных условий), однако для масштабов моделирования (десятки, сотни тысяч частиц) объемное охлаждение реализуется достаточно просто, а результаты таких компьютерных экспериментов хорошо согласуются с эмпирикой.

Более сложным в реализации является проведение моделирования процесса изобарного охлаждения. В этом случае предполагается использование процедур корректировки скоростей частиц и баростатирования. Баростатирование должно выполнять изменение размеров ячейки моделирования и расстояния между частицами для поддержания постоянного давления системы. На практике применяются два основных метода поддержания давления: баростат Андерсена [20] и баростат Берендсена [21]. Метод Андерсена коренным образом меняет алгоритм Верле путем введения обобщенных координат, связанных с размерами системы и скоростью изменения этих размеров. Этот алгоритм довольно сложно реализуем, а жесткость баростата не позволяет моделировать процессы быстрого изменения температуры. Для работы баростата Берендсена необходимо использовать табличные значения изотермической сжимаемости для каждой точки фазовой траектории, что для моделирования неравновесных процессов невозможно.

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

Так, например, теплоемкость Cυ связана с флуктуациями температуры

$\frac{{\left\langle {{{T}^{2}}} \right\rangle - {{{\left\langle T \right\rangle }}^{2}}}}{{{{{\left\langle T \right\rangle }}^{2}}}} = \frac{3}{{2N}}\left( {1 - \frac{{3kN}}{{2{{С }_{\upsilon }}}}} \right),$
а на основе известной парной корреляционной функции распределения h(r) можно определить изотермическую сжимаемость системы

$kT\left( {\frac{{dn}}{{dp}}} \right) = 1 + 4\pi n\int\limits_0^\infty {h\left( r \right){{r}^{2}}dr} .$

Способы расчета калориметрических величин по тепловым флуктуациям резко критикуются в работах В.Ю. Бугаева [22], так как они справедливы в случае, когда усреднение по времени соответствует усреднению по каноническому ансамблю, что для моделирования процессов абсолютно неприемлемо. Так же недостоверны расчеты свойств системы по данным о парной корреляционной функции в условиях протекания интенсивного термодинамического процесса – для получения корректной функции распределения необходимо усреднение по большому числу итераций в равновесном или квазиравновесном состояниях.

В [23] описан один из способов моделирования ансамбля с постоянным числом частиц, давлением и температурой. Здесь используется баростат, корректирующий положения частиц и размеры системы на коэффициент µ, выведенный из лагранжиана (NPT)-ансамбля:

$\mu = 1 + \frac{{P(V) - P}}{{3P(V) + \sum\limits_{i < j} {{{{\text{Ф }}}_{{ij}}}{{r_{{ij}}^{2}} \mathord{\left/ {\vphantom {{r_{{ij}}^{2}} {3V}}} \right. \kern-0em} {3V}}} }},$
где P – заданное давление; P(V) – давление, рассчитанное на данной итерации моделирования. Такой коэффициент пригоден для небольших флуктуаций давления и практически не поддерживает постоянным давление в процессе быстрого охлаждения. Однако этот метод можно применить в циклической процедуре (рис. 1), позволяющей поддерживать давление с точностью 1% даже в случаях резких перепадов температуры.

Рис. 1.

Алгоритм процедуры баростатирования.

Моделирование изобарного охлаждения дает возможность численного расчета таких характеристик, как изобарная теплоемкость Сp и изотермическая сжимаемость βT, динамика которых позволяет прослеживать закономерности фазовых переходов II рода.

Так, получив в результате компьютерного эксперимента температурную зависимость энтальпии H, можно численно определить изобарическую теплоемкость, опираясь на физический смысл этой величины:

${{C}_{p}} = {{\left( {\frac{{dH}}{{dT}}} \right)}_{P}}.$

В разработанной методике расчета изотермической сжимаемости необходимо получить две изобары с близкими значениями давления P1 и P2. По разностям плотностей на фазовых кривых изотермическая сжимаемость определяется следующим выражением:

${{\beta }_{Т }} = \frac{1}{{\left\langle \rho \right\rangle }}\frac{{{{\rho }_{{P1}}} - {{\rho }_{{P2}}}}}{{{{P}_{1}} - {{P}_{2}}}},$
где ρP1 и ρP2 – плотности систем под давлением P1 и P2 при одинаковых температурах.

Моделирование процесса адиабатического расширения дает возможность численного определения адиабатической сжимаемости βS. Предлагаемая методика моделирования такого процесса включает следующие этапы.

1. На протяжении заданного количества итераций моделирования n производится усреднение давления 〈P1〉 и внутренней энергии системы 〈E1〉 (для времени моделирования Δt = 10–15 с количество итераций усреднения – порядка 1000).

2. Размеры ячейки моделирования увеличиваются, как и расстояние между частицами, на величину dV, характеризующую скорость расширения.

3. Определяются средние производной вириала, кинетическая 〈Ek2〉 и потенциальная 〈U2〉 энергии полученной системы для следующих n итераций моделирования.

4. Производится корректировка скоростей системы на коэффициент k, выводимый из нижеследующих соображений.

Запишем условие отсутствия притока–оттока тепла

(1)
$0 = \left\langle {E_{{k2}}^{*}} \right\rangle + \left\langle {{{U}_{2}}} \right\rangle - \left\langle {{{E}_{1}}} \right\rangle - \frac{{{{P}_{1}} + {{P}_{2}}}}{2}dV.$

Здесь $\left\langle {E_{{k2}}^{*}} \right\rangle $ – средняя кинетическая энергия, которой должна обладать система после расширения для выполнения условия постоянства энтропии. В свою очередь, давление P2 определяется из теоремы о вириале

(2)
${{P}_{2}} = \frac{\rho }{3}\left( {2\left\langle {E_{{k2}}^{*}} \right\rangle - \frac{1}{N}\left\langle {\sum\limits_{i = 1}^{N - 1} {\sum\limits_{j > i} {{{r}_{{ij}}}\frac{{\partial \Phi }}{{\partial {{r}_{{ij}}}}}} } } \right\rangle } \right).$

Подставив (2) в (1), выразим $\left\langle {E_{{k2}}^{*}} \right\rangle $ и выполним нормировку ${{\left\langle {E_{{k2}}^{*}} \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {E_{{k2}}^{*}} \right\rangle } {\left\langle {{{E}_{{k2}}}} \right\rangle }}} \right. \kern-0em} {\left\langle {{{E}_{{k2}}}} \right\rangle }}.$ Соответствующий нормировочный коэффициент для скоростей частиц системы определяется выражением

(3)
$k = \sqrt {\frac{{\left\langle {{{E}_{1}}} \right\rangle - \left\langle {{{U}_{2}}} \right\rangle - \frac{{dV}}{2}\left( {\left\langle {{{P}_{1}}} \right\rangle - \frac{\rho }{{3N}}\left\langle {\sum\limits_i {\sum\limits_{j > i} {{{r}_{{ij}}}\frac{{\partial {\text{Ф }}}}{{\partial {{r}_{{ij}}}}}} } } \right\rangle } \right)}}{{\left( {1 + \frac{\rho }{{3N}}} \right)\left\langle {{{E}_{{k2}}}} \right\rangle }}} .$

5. Повторяем пп. 1–4 до установления необходимого объема.

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

Стоит отметить, что данная методика пригодна для моделирования медленного расширения, так как в случае больших значений dV выражение (3) перестает отвечать требованиям корректности численного интегрирования.

По значениям объема и давления, полученным в результате моделирования, определяется адиабатическая сжимаемость

${{\beta }_{S}} = - \frac{1}{V}{{\left( {\frac{{dV}}{{dP}}} \right)}_{S}}.$

Теплоемкости и сжимаемости вещества для фазовой точки связаны между собой соотношением

(4)
${{\beta }_{S}} = \frac{{{{C}_{v}}}}{{{{C}_{p}}}}{{\beta }_{T}}.$

Сведения об этих величинах позволяют описать всю термодинамику системы, определить значения упругих характеристик (объемный модуль упругости, скорость звука и др.), а при известной парной корреляционной функции получить термическое уравнение состояния для моделируемой системы частиц.

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

На основе описанных выше методик разработана программа для моделирования молекулярных систем MDDX11, на которую получено свидетельство о государственной регистрации [24]. С помощью этой программы проведен ряд экспериментов, моделирующих процессы равновесного и неравновесного изобарического охлаждения и адиабатического расширения систем частиц аргона.

На рис. 2 представлены результаты расчета сжимаемости по данным моделирования равновесного и неравновесного охлаждения соответственно.

Рис. 2.

Изотермическая сжимаемость аргона при 4 МПа: (а) 1 – результаты компьютерных экспериментов, моделирующих охлаждение аргона со скоростью 109 К/с; 2 – данные [25]; (б) – результаты моделирования со скоростью охлаждения 1012 К/с.

На графике температурной зависимости сжимаемости (рис. 2а), полученной по данным эксперимента с равновесным медленным охлаждением, можно различить газовую, жидкую и кристаллическую фазы.

На рис. 2б представлен ряд полученных значений β при охлаждении аргона со скоростью 1012 К/с. Температурная зависимость сжимаемости β испытывает резкое изменение в окрестностях ~90 К. Такое резкое изменение, видимо, характеризует неравновесный процесс при быстром охлаждении моделируемой системы частиц аргона. Возможно, это связано с переходом “неравновесный газ (перенасыщенный пар)–неравновесная жидкость (смесь газа с кластерами жидкости)”. Также можно отметить переход в окрестностях ~30 К, вероятно соответствующий переходу нестабильной жидкости в стеклообразную фазу.

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

Рис. 3.

Адиабатическая сжимаемость аргона при 4 МПа по данным компьютерных экспериментов, моделирующих адиабатическое расширение в окрестностях фазовых точек изобары.

Полученные значения адиабатической сжимаемости отвечают зависимости параметров, описываемой выражением (4).

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

ЗАКЛЮЧЕНИЕ

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

Работа выполнена при поддержке РФФИ (грант № 18-42-030002 р_а).

Список литературы

  1. Сандитов Б.Д., Сангадиев С.Ш., Цыдыпов Ш.Б. Взаимодействие между параметром Грюнайзена и коэффициентом Пуассона стеклообразных систем // Акуст. журн. 2007. Т. 53. № 4. С. 112.

  2. Бадмаев Б.Б., Бальжинов С.А., Дамдинов Б.Б., Дембелова Т.С. Низкочастотная сдвиговая упругость жидкостей // Акуст. журн. 2010. Т. 56. № 5. С. 602.

  3. Есипов И.Б., Зозуля О.М., Миронов М.А. Медленная кинетика нелинейности вязкоупругих свойств нефти при сдвиговых колебаниях // Акуст. журн. 2014. Т. 60. № 2. С. 166.

  4. Мартынов Г.А. Классическая статистическая механика. Теория жидкостей. Долгопрудный: Интеллект, 2011. 328 с.

  5. Trozzi G., Ciccoti G. Stationary Nonequilibrium States by Molecular Dynamics. II. Newton’s Law // Phys. Rev. 1984. V. 29. № 1. P. 916.

  6. Evans D.J. Thermal Conductivity of the Lennard-Jones Fluid // Phys. Rev. 1984. V. 29. № 1. P. 1449.

  7. Paolini G.V., Ciccoti G., Massobrio C. Nonlinear Thermal Response of a Lennard-Jones Fluid near the Triple Point // Phys. Rev. A – Gen. Phys. 1986. V. 34. № 2. P. 1355.

  8. Hoover W.G. Computer Simulation of Many-body Dynamics // Phys. Today. 1984. V. 37. № 1. P. 44.

  9. Polymeropoulos E.E., Brickman J. Molecular Dynamics Study of the Formation of Argon Clusters in the Compressed Gas // Chem. Phys. Lett. 1982. V. 92. № 1. P. 59.

  10. Желудков С.В., Инсепов 3.А. Молекулярно-динамическое моделирование кинетики конденсации в пересыщенном паре // ЖФХ. 1987. Т. 61. № 4. С. 1109.

  11. Yonezawa F., Nose Sh., Sakamoto Ch. Computer Simulation of the Glass Transition // Proc. VI Int. Conf. on Liquid and Amorphous Metals. V. 1. Garmish-Partenkirchen. 1986. Oldenburg, 1987. P. 77.

  12. Zarzychi J. Glass Structure // J. Non-Cryst. Solids. 1982. V. 52. № 1/3. P. 31.

  13. Nagel S.R., Grest G.S., Rahman A. Quench Echoes // Phys. Today. 1983. V. 36. № 10. P. 24.

  14. Angell C.A., Clarke J.H.R., Woodcock L.V. Interaction Potentials and Glass Formation: a Survey of Computer Experiments // Adv. Chem. Phys. 1981. V. 48. P. 397.

  15. Jackle J. Models of the Glass Transition // Rep. Prog. Phys. 1986. V. 49. P. 171.

  16. Kimura M., Yonezava F. Computer Glass Transition // Topological Disorder in Condensed Matter: Proc. of the 5 Taniguchi Int. Symp. / Ed. Yonezava F., Ninomiya T. Berlin: Springer Verlag, 1983. P. 80.

  17. Цыдыпов Ш.Б., Парфенов А.Н., Сандитов Д.С., Аграфонов Ю.В., Нестеров А.С. Применение метода молекулярной динамики и модели возбужденного состояния к изучению процесса стеклования аргона // Физ. и хим. стекла. 2006. Т. 32. № 1. С. 116.

  18. Цыдыпов Ш.Б., Сандитов Д.С. Критерий стеклования жидкостей в модели возбужденных атомов // ЖФХ. 2004. Т. 78. № 5. С. 906.

  19. Цыдыпов Ш.Б., Сандитов Д.С., Парфенов А.Н. Исследование стеклования аргона методом молекулярной динамики // ЖФХ. 2005. Т. 79. № 9. С. 1464.

  20. Andersen H.C. Molecular Dynamics Simulations at Constant Pressure and/or Temperature // J. Chem. Phys. 1980. V. 72. P. 2384.

  21. Berendsen H.J.C. Molecular Dynamics with Coupling to an External Bath // J. Chem. Phys. 1984. V. 81. P. 3684.

  22. Бугаев В.Ю., Рабинович В.А. О методах расчета термодинамических свойств жидкостей в молекулярно-динамических экспериментах // ТВТ. 1983. Т. 5. № 21. С. 871.

  23. Rapaport D.C. The Art Molecular Dynamics Simulation. Cambridge: Cambridge Univ. Press, cop., 2004. 605 p.

  24. Герман Е.И. Программа моделирования молекулярных систем MDDX11. Св-во о гос. рег. программы для ЭВМ № 2 016 617 783.

  25. Александров А.А., Орлов К.А., Очков В.Ф. Свойства и процессы рабочих тел и материалов атомной энергетики. Интеракт. спр. http://twt.mpei.ac.ru/TTHB/ NPP/

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Теплофизика высоких температур

Архивы выпусков Информация о журнале Отправить рукопись в журнал