CC BY
76
УДК 538.931:544.272.23:546.72
ВЯЗКОСТЬ ЖИДКОГО ЖЕЛЕЗА: МОЛЕКУЛЯРНО-ДИНАМИЧЕСКИЙ РАСЧЕТ С ПОТЕНЦИАЛОМ ПОГРУЖЕННОГО АТОМА
И.В. Мальцев, А.А. Мирзоев
Метод молекулярной динамики с соотношением Грина-Кубо применен для расчета вязкости жидкого железа в широком интервале температур вблизи температуры плавления. Показано существование особенностей в температурной зависимости вязкости. Сделан вывод о связи внутренней структуры модели и поведения вязкости.
Ключевые слова: молекулярная динамика, сдвиговая вязкость, метод погруженного атома, жидкое железо.
Введение
Железо - доминирующий компонент большинства используемых человечеством металлических конструкций. Механические свойства железа и его сплавов определяются предысторией расплава перед кристаллизацией и условиями кристаллизации [1]. Были обнаружены гистерезис свойств (вязкость, поверхностное натяжение, магнитная восприимчивость, плотность) при нагреве и охлаждении расплава, осцилляция этих свойств со временем. В настоящее время отсутствует общепринятое достоверное описание связи свойств расплава и его внутренней структуры. Экспериментальное исследование структуры и свойств расплавов на основе железа затруднено как в силу ряда технических причин, так и из-за чувствительности расплавов к концентрации компонентов. Молекулярно-динамическое (МД) моделирование [2] с использованием потенциалов погруженного атома [3] позволяет строить модели, соответствующие реальным системам, в частности, в работе [4] был разработан ЕАМ (Embedded atom method, метод погруженного атома) потенциал для жидкого железа, хорошо предсказывающий температуру плавления, плотность, энергию. В работах [5, 6] предложены различные версии потенциала для железа при больших давлениях и температурах и применены соответственно в [7] и [6] для расчета вязкости в условиях, близких к условиям в ядре Земли.
Данные о температурной зависимости вязкости железа вблизи температуры плавления отсутствуют в настоящее время. Представляется интересным восполнение этого пробела и соотнесение молекулярно-динамических расчетов с экспериментальными данными. Наличие или отсутствие особенностей на функции вязкости может свидетельствовать за или против выбранного взаимодействия, представленного ЕАМ потенциалом. Кроме того, анализ и интерпретация результатов могут приблизить нас к пониманию внутренней структуры расплава железа.
В работе представлен расчет сдвиговой вязкости чистого железа с ЕАМ потенциалом [4]. Расчет проводился в рамках метода равновесной молекулярной динамики с использованием метода Грина-Кубо [8, 9].
Методика расчета
Кинетические коэффициенты, в частности коэффициент сдвиговой вязкости, могут быть определены в рамках теории линейного отклика методом Грина-Кубо [8, 9]. Коэффициент переноса выражается как интеграл от автокорреляционной функции соответствующей величины (1). Известно, что для вычисления коэффициента сдвиговой вязкости необходимо найти автокорреляционную функцию тензора давления:
¥
h = j x(t)dt, (1)
о
где h - коэффициент сдвиговой вязкости, X (t) - автокорреляционная функция тензора давления (2),
x(t)=(0У Pa (tУ • (2)
где V - объем системы, T - температура системы, Pap (t) - недиагональная компонента тензора давления (3),
Физика
N
Іар() = V X miVia і=і
N
і< і
(3)
где N - число частиц, тг- - масса / -й частицы, \1а (7) - а -я компонента скорости / -й частицы, Г1]а () - а -я компонента расстояния между частицами /', /, ¥^р () - сила, действующая между частицами /', /.
В формуле (2) проводится усреднение по ансамблю. Для получения качественного усреднения необходимо проводить расчет для большого количества (> 10 000) различных исходных конфигураций. Подготовка и расчет такого набора независимых конфигураций занимает существенное вычислительное время, поэтому для оптимизации временных затрат была выбрана следующая схема. Проводится расчет одной исходной конфигурации с количеством шагов интегрирования, равным 10 000 000 (шаг интегрирования уравнений движения -
0,001 пс). Далее из всей полученной траектории через каждые 100 шагов выбираются субтраектории длиной
1 000 000 шагов (рис. 1). Каждая из них считается независимой и участвует в построении автокорреляционной функции. Таким образом, получаем 90 000 траекторий для получения среднего по ансамблю. Для улучшения качества автокорреляционной функции выполнялся расчет от 40 до 160 исходных конфигураций.
Расчеты, в которых непосредственно получался тензор давления, проводились в КУТ ансамбле. Среднее давление было близко к нулю. К кубической ячейке с 2000 атомов железа, взаимодействующих посредством ЕАМ потенциала [4], были приложены периодические граничные условия. Исходные конфигурации были получены двумя способами.
Первый способ
1. Построение одной конфигурации из 250 частиц для твердой фазы.
2. Нагревание до 2500 К полученной конфигурации.
3. Постепенное охлаждение с запоминанием конфигураций при требуемых для исследования температурах.
4. Репликация конфигураций в каждом направлении дважды (2х2х2); таким образом, получаем по одной конфигурации из 2000 частиц для заданной температуры.
5. Для каждой температуры многократно пересоздаем скорости (из нормального распре -деления) и проводим уравновешивание; таким образом, получаем требуемое количество конфигураций для каждой температуры.
Второй способ
1. Построение требуемого количества конфигураций из 2000 частиц для твердой фазы.
2. Нагревание до 2200 К каждой конфигурации.
3. Постепенное охлаждение с запоминанием конфигураций при требуемых для исследования температурах.
Рис. 1. Схема построения автокорреляционной функции
Результаты моделирования
Молекулярно-динамический расчет проводился для температур из интервала 1650-2100 К. Температура плавления чистого железа составляет около 1812 К, следовательно, интервал 16501812 К соответствует переохлажденной жидкости. Для верификации расчета была получена плотность жидкого железа в МД модели. Результат с известными экспериментальными данными изображен на рис. 2. При температурах, выше точки плавления, МД значения плотности пре-
имущественно лежат в пределах области разброса экспериментальных данных [10, с. 144]. Виден
различный наклон функций плотности. Хорошее совпадение численных значений плотности свидетельствует о пригодности выбранного потенциала [4] и отсутствии ошибок в МД расчете.
Коэффициент сдвиговой вязкости (рис. 3) рассчитан для двух наборов конфигураций, как описано в методике расчета. Как для первого набора, так и для второго есть области плавного и скачкообразного изменения вязкости. Интервал температур от 1766 до 1816 К соответствует переохлажденной жидкости. Значения вязкости здесь изменяются немонотонно, причем очень высок разброс значений. Это сигнали-Рис. 2. Температурная зависимость плотности жидкого железа: зирует о возможной нестабильности
квадратные маркеры - эксперимент [10], круглые - МД расчет „
внутренней структуры расплава в
этой области. Значения вязкости для первого набора получены с большим усреднением (160 конфигураций), тем не менее провал на функции вязкости остается хорошо заметным. Ниже температуры 1750 К вязкость изменяется плавно, что скорее всего указывает на плавное изменение параметров структуры. Выше температуры плавления в каждом из наборов конфигураций наблюдается участок (1833-1883 К - набор 1, 1816-1850 К - набор 2) с плавным изменением вязкости. Далее (до 1950 К) в каждой из зависимостей наблюдаются скачкообразные изменения вязкости, а также почти горизонтальные участки. Это сигнализирует об интенсивном изменении
,, ------данные 12 [10, стр. 146]
__I__I_I__I_I__I_I__I_I__I__I_I__I_I__I I__I_I__■ I__I__I_I__I_I__I_I_____I_I__I_I_^
1700 1750 1800 1850 1900 1950 2000
Т,К
Рис. 3. Температурная зависимость вязкости: круглые маркеры - набор конфигураций 1; треугольные - набор 2; линии без маркеров - экспериментальные данные (см. легенду)
Физика
структуры в этом интервале температур, что отмечалось во множестве экспериментальных работ [1]. Одинаковое поведение вязкости в этом интервале для рассмотренных наборов указывает на
существование определенной структуры, обусловливающей обнаруженный характер температурной зависимости. Простое косвенное подтверждение изменения структуры можно найти в изменении высоты второго пика парной корреляционной функции с температурой (рис. 4). В указанных выше интервалах наблюдается немонотонное или резкое изменение высоты пика в сторону увеличения, что говорит о существовании упорядочения внутренней структуры. Это также подтверждает связь внутренней структуры и поведения вязкости.
Основное различие исследованных наборов конфигураций заключается в смещении интервалов с особенностями. Такое различие объясняется отличающимися условиями формирования исходных моделей, однако явная связь между условиями образования конфигураций и видом вязкости пока остается неясной.
Также на рис. 3 показаны некоторые экспериментальные результаты. Полученные значения вязкости лежат между экспериментальными кривыми, что свидетельствует о пригодности используемого потенциала [4] для расчета вязкости жидкого железа вблизи точки плавления.
Выводы
В этой работе применен метод молекулярной динамики с соотношением Грина-Кубо для изучения температурной зависимости вязкости жидкого железа. Для построения модели был использован потенциал погруженного атома, разработанный Менделевым и др. [4]. Показано, что этот потенциал корректно воспроизводит как плотность, так и вязкость жидкого железа. Обнаружены особенности зависимости вязкости от температуры: скачкообразные изменения и области медленного изменения вязкости (горизонтальные участки). Наиболее вероятной причиной этих особенностей предполагается наличие упорядочения в структуре и переход от одного упорядочения к другому. Более строгое обоснование этого предположения должно быть связано с исследованием структуры моделей и соотнесением её с характером изменения вязкости.
Работа поддержана грантом РФФИ № 09-03-00584.
Литература
1. Свойства металлических расплавов: сб. научн. тр.: в 2 ч. / В.С. Цепелев, В.В. Конашков, Б.А. Баум и др. - Екатеринбург: УГТУ-УПИ, 2008. - Ч. 1.- 358 с.
2. Allen, M.P. Computer Simulation of Liquids / M.P. Allen, D.J. Tildesley. - Oxford Univ. Press.: Oxford, 1987. - 385 c.
3. Daw, M.S. Semiempirical, Quantum Mechanical Calculation of Hydrogen Embrittlement in Metals / M.S. Daw, M.I. Baskes // Phys.Rev.Lett. - 1983. - V. 50. - P. 1285-1288.
4. Mendelev, M.I. / M.I. Mendelev, S. Han, D.J. Srolovitz et al.// Phil.Mag.A. - 2003. - V. 83. -P. 3977.
5. Koci, L. / L. Koci, A.B. Belonoshko, R. Ahuja // Phys.Rev.B. - 2006. - P. 224113.
6. Белащенко, Д.К. Применение модели погруженного атома к жидким металлам. Жидкое железо / Д.К. Белащенко // Журнал физической химии. - 2006. - V. 80. - P. 1-11.
7. Desgranges, C. Viscosity of liquid iron under high pressure and high temperature: Equilibrium and nonequilibrium molecular dynamics simulation studies / C. Desgranges, J. Delhommelle //Phys.Rev.B. - 2007. - V. 76. - P. 172102.
8. 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. 398413.
9. 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.
10. Еланский Г.Н. Строение и свойства металлических расплавов / Г.Н. Еланский, Д.Г. Еланский. - М.: МГВМИ, 2006. - 228 с.
Поступила в редакцию 5 июня 2009 г
LIQUID IRON VISCOSITY: MOLECULAR-DYNAMICS SIMULATION WITH AN EMBEDDED-ATOM POTENTIAL
The shear viscosity of liquid iron was calculated through equilibrium (Green-Kubo) molecular-dynamics simulation in a wide range of temperatures near the melting point. It was found that there are specific features in the temperature dependence of the shear viscosity. It was pointed out that there is a structure-viscosity relationship.
Keywords: molecular dynamics, shear viscosity, embedded-atom method, liquid iron.
Maltsev Ilya Vladimirovich - Post-Graduate Student, General and Theoretical Physics Department, South Ural State University.
Мальцев Илья Владимирович - аспирант, кафедра общей и теоретической физики, ЮжноУральский государственный университет.
e-mail: maltsev.ilya@gmail.com
Mirzoev Aleksandr Aminulaevich - Dr.Sc. (Physics and Mathematics), Professor, General and Theoretical Physics Department, South Ural State University.
Мирзоев Александр Аминулаевич - доктор физико-математических наук, профессор, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
e-mail: mirzoev@physics.susu.ac.ru
