CC BY
53
УДК 538.931:544.272.23:546.72
МЕТОДИКА РАСЧЕТА СДВИГОВОЙ ВЯЗКОСТИ ЖИДКОГО ЖЕЛЕЗА МЕТОДОМ ГРИНА-КУБО
И.В. Мальцев, А.А. Мирзоев
Рассмотрен вопрос о применении равновесной молекулярной динамики к расчету коэффициента сдвиговой вязкости жидкого железа. Для расчета использован метод Грина-Кубо. Построены зависимости вязкости от параметров используемого термостата, количества частиц системы, шага интегрирования.
Ключевые слова: молекулярная динамика, сдвиговая вязкость, жидкое железо, метод Грина-Кубо.
Введение
Метод молекулярной динамики уже давно используется для исследования разнообразных физических систем. В последнее десятилетие были разработаны [1-3] потенциалы взаимодействия, позволяющие моделировать системы атомов железа в различных условиях. Один из наиболее подходящих потенциалов для моделирования жидкого железа при нормальных условиях был предложен Менделевым и др. [2].
Ранее уже рассматривались вопросы о влиянии параметров моделирования и точности вычисления коэффициента сдвиговой вязкости [4, 5] для жидкости Леннарда-Джонса. Было показано [4], что коэффициент сдвиговой вязкости (далее вязкость) очень слабо зависит от числа частиц. Для количества частиц больше 256 вычисленные значения лежат близко друг к другу в пределах погрешности. Верле и др. [6] не обнаружили зависимости от числа частиц.
В литературе практически отсутствует информация по влиянию внутренних параметров моделирования: термостат, шаг интегрирования. Для проведения масштабных численных экспериментов по определению вязкости жидкого железа представляется важным рассмотрение этих вопросов и нахождение параметров, при которых вышеупомянутая зависимость от числа частиц отсутствует.
Метод расчета
Вязкость в настоящей работе рассчитывалась методом Грина-Кубо [7, 8], для чего приведенная в равновесие система частиц, взаимодействующих посредством потенциала [4], развивалась продолжительной время. Подробности применения метода Грина-Кубо могут быть найдены в работе [9]. Температура в системе была установлена равной T = 1950 K, давление равнялось нулю. Широко известный термостат Нозье-Гувера [10], использованный в этой работе, имеет один важный варьируемый параметр, который отвечает за скорость релаксации температуры. Пакет LAMMPS [11], при помощи которого проводилось моделирование, принимает этот параметр как время, необходимое для релаксации системы.
В табл. 1 приведены параметры проведенных численных экспериментов. Каждый эксперимент включал в себя 100 независимых расчетов с различными начальными конфигурациями. В каждом эксперименте для получения автокорреляционной функции тензора давления участвовало 39 000 субтраекторий [9], и только при шаге интегрирования 0,0001 10-12 усреднение проводилось по 18 500 субтраекториям.
Коэффициент самодиффузии был рассчитан методом Эйнштейна-Гельфанда [12]. Так как вязкость в методе Грина-Кубо представлена как интеграл по времени от автокорреляционной функции тензора давления, то необходимо задать время, при котором следует оборвать вычисление интеграла. Мы использовали здесь значение времени, достаточное, чтобы интеграл изменялся мало (около некоторого среднего), что автоматически приводит к учету медленно сходящейся к нулю части автокорреляционной функции при больших t («long-time tail»).
Это время было принято равным 10-11 с, что соответствует 10 000 шагов моделирования при шаге, равном 0,001 10-12 с. Определение погрешности An(t) было основано на хорошо известной формуле для стандартного отклонения от среднего п* (t), примененной в работе [5]:
Мальцев И.В., Мирзоев А.А.
Методика расчета сдвиговой вязкости жидкого железа методом Грина-Кубо
г ) =
(г)-п (г))2
_к_
N (N -1)
'ехр^ 'ехр у
12
(1)
где Nexp - количество независимых экспериментов, пк(г) - вязкость, полученная в отдельном эксперименте.
Таблица 1
Количество частиц N Шаг интегрирования Л,10-12 с Параметр термостата гге1ах ,10 12 С Коэффициент самодиффузии Д10-9 м2/ с Коэффициент сдвиговой вязкости -2 П,10 пуаз
2000 0,001 10,0 4,28 ± 0,02 4,15 ± 0,02
4000 0,001 10,0 4,38 ± 0,02 4,17 ± 0,01
2000 0,0001 10,0 4,27 ± 0,03 4,17 ± 0,03
2000 0,0002 10,0 4,27 ± 0,05 4,15 ± 0,02
2000 0,001 1,0 4,29 ± 0,02 4,13 ± 0,02
2000 0,001 100,0 4,27 ± 0,02 4,18 ± 0,02
2000 0,001 1000,0 4,27 ± 0,03 4,14 ± 0,02
Результаты и обсуждение
Из табл. 1 видно, что коэффициент самодиффузии монотонно уменьшается при увеличении времени релаксации для термостата, но при этом разница между вычисленными значениями менее 1 %. Коэффициент сдвиговой вязкости не обнаруживает (рис. 1) монотонной зависимости от параметра термостата, разница значений здесь также мала и составляет около 1 %. Этот результат позволяет при выборе параметра термостата руководствоваться минимизацией времени для приведения и удержания модели в равновесии.
Значения вязкости систем с различным числом частиц также очень близки, что подтверждает для моделируемой системы жидкого железа результат, полученный для жидкости Леннарда-Джонса. Так как компьютерный расчет вязкости занимает продолжительное время, то возможность ограничиваться малым (~2000 частиц) размером системы позволяет обеспечить лучшее усреднение автокорреляционной функции при той же длительности моделирования.
Отсутствие практически значимой разницы коэффициента вязкости при различном выборе шага интегрирования также позволяет использовать большой шаг интегрирования. Для оценки
верхней границы шага интегрирования, при котором еще возможно получение коэффициента вязкости, был проведен молекулярно-динамический эксперимент в шагом
0,005 10-12 с. Было обнаружено, что система частиц железа перестает быть устойчивой: расчет вязкости дает заведомо неверный результат.
Получение корректной автокорреляционной функции тензора давления требует усреднения по большому количеству начальных конфигураций. От ее точности во многом зависит сходимость интеграла вязкости в методе Грина-Кубо. Использование в работе 100 независимых конфигураций позволило добиться погрешности расчета вязкости менее 1 %.
Полученное в настоящей работе значение вязкости совпадает со значением, полученным ранее [9], что подтверждает воспроизводимость расчетов.
4.2
4.19
4.18
4.17
4.16
* 4.15
4.13
4.11
О
о
о
о
10
10
10
10
Рис. 1. Зависимость коэффициента сдвиговой вязкости от параметра термостата
Серия «Математика. Механика. Физика», выпуск 2
77
Физика
Заключение
Проведенные численные эксперименты показали слабую зависимость коэффициента сдвиговой вязкости от параметров моделирования: разница значений составляет около 1%. Выявлен предельный шаг интегрирования 0,005 10-12 с, при котором система частиц становится неустойчивой и не дает значения вязкости. Полученные данные могут быть применены для планирования молекулярно-динамических расчетов вязкости.
Численные эксперименты проведены на базе суперкомпьютерного центра ЮУрГУ [13].
Литература
1. Mendelev, M.I. Development of New Interatomic Potentials Appropriate for Crystalline and Liquid Iron / M.I. Mendelev, S. Han, D.J. Srolovitz et al. // Phil.Mag. A. - 2003. - V. 83. - P. 3977.
2. Koci, L. Molecular dynamics study of liquid iron under high pressure and high temperature / L. Koci, A.B. Belonoshko, R. Ahuja // Phys. Rev. B. - 2006. - V. 73. - P. 224113.
3. Белащенко, Д.К. Применение модели погруженного атома к жидким металлам. Жидкое железо / Д.К. Белащенко // Журнал физической химии. - 2006. - V. 80. - P. 1-11.
4. Jerome J. Erpenbeck Shear viscosity of the Lennard-Jones fluid near the triple point: Green-Kubo results / Erpenbeck Jerome J. // Phys. Rev. A. - 1988. - V. 38. - P. 6255.
5. Ryckaert, J.-P. Evaluation of transport coefficients of simple fluids by molecular dynamics: om-parison of Green-Kubo and nonequilibrium approaches for shear viscosity / J.-P. Ryckaert, A. Belle-mans, G. Ciccotti, G.V. Paolini // Phys. Rev. A. - 1989. - V. 39. - P. 259.
6. Levesque, D. Computer «Experiments» on Classical Fluids. IV. Transport Properties and Time-Correlation Functions of the Lennard-Jones Liquid near Its Triple Point / D. Levesque, L. Verlet, J. Kurkijarvi // Phys. Rev. A. - 1973. - V. 7. - P. 1690.
7. Green, M.S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids / M.S. Green // J. Chem. Phys.-1954.-V.22.-P.398-413.
8. Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems / R. Kubo // J. Phys. Soc. Jpn. - 1957. - V. 12. -P. 570-586.
9. Мальцев, И.В. Вязкость жидкого железа: молекулярно-динамический расчет с потенциалом погруженного атома / И.В. Мальцев, А.А. Мирзоев // Вестник ЮУрГУ. Серия «Математика. Механика. Физика». - 2009. - Вып. 1. - № 22(155). - С. 79-83.
10. Nose, S. A unified formulation of the constant temperature molecular dynamics methods / S. Nose//J. Chem. Phys. - 1984. -V. 81. - P. 511.
11. http://lammps.sandia.gov/
12. Helphand, E. Transport Coefficients from Dissipation in a Canonical Ensemble / E. Helphand // Phys. Rev. - 1960. -V. 119. - P. 1.
13. http://supercomputer. susu.ru/computers/ckif_ural/
Поступила в редакцию 8 февраля 2010 г.
APPLICATION OF GREEN-KUBO METHOD TO LIQUID IRON SHEAR VISCOSITY CALCULATION
Green-Kubo method within equilibrium molecular-dynamics simulation was applied to calculate shear viscosity of liquid iron. Simulation-parameters (thermostat strength, number of particles, integration step) dependence of shear viscosity coefficient was obtained.
Keywords: molecular-dynamics simulation, shear viscosity, liquid iron, Green-Kubo method.
Maltsev Ilya Vladimirovich - post-graduate student, General and Theoretical Physics Department, South Ural State University.
Мальцев Илья Владимирович - аспирант, кафедра общей и теоретической физики, ЮжноУральский государственный университет.
e-mail: maltsev.ilya@gmail.com.
Mirzoev Alexander Aminulaevich - Dr.Sc. (Physics and Mathematics), Professor, General and Theoretical Physics Department, South Ural State University.
Мирзоев Александр Аминулаевич - профессор, доктор физико-математических наук, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
