CC BY
108
УДК 548.3
Г.М. Полетаев, А.М. Сагалаков, П. С. Стенченко Молекулярно-динамическое моделирование и стохастические свойства аморфных металлов
G.M. Poletaev, A.M. Sagalakov, P.S. Stenchenko
Molecular Dynamics Simulation and Stochastic Properties of Amorphous Metals
С позиции теории динамических систем рассматривается вопрос моделирования методом молекулярной динамики аморфных металлов. Показано, что для молекулярно-динамических моделей характерны неустойчивость по Ляпунову и структурная неустойчивость по Понтрягину, в связи с чем конкретная атомная структура аморфных металлов априори является непредсказуемой и невоспроизводимой.
Ключевые слова.: молекулярная динамика, аморфный металл, неустойчивость по Ляпунову, неустойчивость по Понтрягину, теория динамических систем, кристаллический кластер.
The question of simulation of amorphous metals by the method of molecular dynamics is considered from the position of the theory of dynamic systems. It is shown that for molecular-dynamics systems is typical the Lyapunov instability and structural instability by Pontryagin. In this connection the concrete atomic structure of amorphous metals is a priory unpredictable and non-reproducible.
Key words: molecular dynamics, amorphous metals, Lyapunov instability, Pontryagin instability, theory of dynamic systems, crystal cluster.
Атомная структура жидких металлов, аморфных металлов и сплавов к настоящему времени сравнительно слабо изучена. Основная трудность - ограниченные возможности прямых экспериментальных методов, которые обусловлены специфическими свойствами изучаемых сред. В то же время аморфные вещества находят все большее применение в различных отраслях науки и техники и уже начинают кое-где вытеснять обычные металлы и сплавы с кристаллической структурой благодаря своим уникальным свойствам: высокой прочности, коррозийной стойкости, изотропности, узкой петли гистерезиса и др. Жидкие металлы используются в качестве теплоносителей на атомных электростанциях, в МГД насосах и генераторах.
Аморфные металлы обладают ярко выраженными металлическими свойствами, но более близки по своей структуре к жидким металлам, из которых они получены путем переохлаждения [1, 2]. Обычно считается, что при очень быстром охлаждении резко увеличивается вязкость жидкости и ее частицы - атомы, молекулы и состоящие из них кластеры - утрачивают подвижность, которую они имели в жидком состоянии. При этом структура жидкости как бы замораживается. Вещество становится твердым в аморфном состоянии, наследуя атомную структуру расплава [1, 2]. Согласно экспериментальным данным, переход от расплава к аморфному твердому телу происходит в определенном диапазоне температур - около температуры стеклования, в котором вязкость расплава очень быстро возрас-
тает, движения атомов затруднены, а характерное время перестроек атомных структур в расплаве сравнивается по порядку величины со временем эксперимента [1].
Однако такая простая качественная картина аморфизации является только некоторым приближением к созданию теории аморфного состояния вещества. Сравнительно легко аморфизируются сплавы сложных систем [1], поэтому изучение атомной структуры расплавов и характера межатомных связей в них играет важную роль в создании теории аморфного состояния.
Структура и свойства полученных металлических стекол зависят от скорости охлаждения [3, 4].
Для изучения атомной структуры аморфных металлов целесообразно использовать метод молекулярной динамики, который состоит в прямом вычислении траекторий атомов в некотором блоке, например, жидкости или твердого тела на основе уравнений движения Ньютона для системы большого числа частиц N с заданным законом межчастич-ного взаимодействия, полученным на основе квантово-механических представлений и содержащим эмпирические постоянные. Далее для получения измеряемых макроскопических величин необходима статистическая обработка - усреднение определенных функций от микроскопических переменных (координат и импульсов частиц системы). Если жидкость или твердое тело подвергаются определенному воздействию, например, нагреваются или охлаждаются, то траектории атомов изменяются,
что приводит к изменению макроскопических свойств системы. Метод молекулярной динамики позволяет получить достаточно полную информацию о физической системе вплоть до частот колебаний атомов и тем самым определить любое свойство системы и даже рассчитывать новые материалы на молекулярном уровне. Данный метод оправдан, если длина волны де Бройля атомов много меньше характерных межатомных расстояний, а температура системы не слишком мала. При очень низких температурах могут стать определяющими квантовые эффекты, и поэтому в данных условиях, вообще говоря, необходимы квантовомеханические расчеты. Традиционно проблемным является вопрос о потенциалах межчастичного взаимодействия, а также об использовании того или иного вида «обрывания» потенциалов.
В данной работе рассматриваются особенности компьютерного моделирования аморфных металлов методом молекулярной динамики с общих позиций теории динамических систем, на которые нередко не обращают никакого внимания (см., однако, работы [5, 6], посвященные обоснованию применимости метода молекулярной динамики для расчетов траекторий отдельных атомов и макроскопических параметров).
В простейшем случае одноатомной среды решают следующую автономную систему дифференциальных уравнений:
dx(m) _ dp(m) ди
dt mt dt dx(m) ’
m _ 1,2,3; i _ 1,..., N; 0 < t < t . (1)
’ ’ ’ ’ ’ > max V '
Здесь N - число атомов в расчетном блоке;
U _ — YUV - потенциальная энергия системы;
2 i * j
U ^ - потенциал взаимодействия i-го и j-го атома.
Это система 6N уравнений, которую необходимо проинтегрировать с начальными условиями
Xi(m>(0) _ xm\ pim) (0) _ p(om). (2)
Если рассматривается процесс охлаждения металла, то задача Коши (1), (2) перестает быть автономной, так как правые части в (1) будут явно зависеть от времени. В этом случае вместо (1), (2) можно использовать динамическую систему с дискретным временем (каскад)
R(n+1) _ Ln)M(n)R(n), i_ 1, 2,..., 6N, где R - б^мерный вектор, включающий в себя компоненты координат и импульсов всех атомов; м(n) - дифференциальный оператор, позволяющий проинтегрировать (1) на одном временном шаге.
Оператор Lf моделирует уменьшение скорости атомов на данном временном шаге за счет охлаждения. Можно использовать, например, линейное убывание скоростей атомов. Неавтономность системы уравнений (1) при охлаждении можно интерпретировать просто как изменение начальных усло-
вий (уменьшение скорости атомов) после каждого временного шага. При такой интерпретации рассматривается совокупность большого числа задач Коши (1), (2) с изменяющимися после каждого временного шага начальными условиями. Отметим, что в рамках данной модели легко рассмотреть и случаи охлаждения с переменной скоростью.
Вернемся к автономной замкнутой системе, моделирующей, например, расплав, в которой полная энергия сохраняется. Давно известно, что динамические системы (1), (2) обычно обладают свойством глобальной устойчивости (финитность движения), которая сочетается с сильной локальной неустойчивостью (неустойчивость движения по Ляпунову) [7]. Задача Коши (1), (2) решается численно, что приводит к неизбежным погрешностям в расчете фазовых траекторий как отдельных атомов, так и всей динамической системы в целом в б^мерном фазовом пространстве. Малые возмущения в процессе численного решения (1) приводят к быстрому отклонению (по экспоненциальному закону) рассчитанной траектории от точного решения при заданных начальных условиях (1). Такое поведение численного решения может быть в принципе проанализировано с помощью первого метода Ляпунова. Расходимость траекторий в данной точке определяется показателем экспоненты - максимальным показателем Ляпунова ^. Усредненная по фазовому пространству величина \ называется энтропией Крылова-Колмогорова - К. Величина К"1 - характерное время данной динамической системы - определяет горизонт предсказуемости будущего системы и одновременно время расцепления корреляций, время релаксации. Величина K в системе одинаковых атомов может быть определена, например, следующим образом [6].
Рассмотрим произвольную начальную точку в фазовом пространстве и будем численно определять выходящие из этой точки фазовые траектории с малым шагом по времени At и другим малым шагом At', например, At' = At /4. Определим величину
(ap2 (t))_ N X (p j (()- pj(t ))2
в совпадающие моменты времени t. Через некоторое время будет выполняться экспоненциальная зависимость
^Ap2 (t)^ _ A exp (Kt) ,
где A -некоторая постоянная.
Далее на основании соотношения
ln ^Ap2 (t )^_ ln A + Kt
определяется величина К. При возрастании времени экспоненциальная зависимость должна смениться насыщением
^Ap2 (t)^ _ 6kTm .
При численном решении по прошествии времени по порядку величины большего 1/К система ато-
мов полностью «забывает» начальные условия. При этом рассчитываемая траектория полностью утрачивает корреляцию с исходной траекторией. Отметим, что, несмотря на указанный выше алгоритм определения К, содержащий свойства вычислительной схемы, энтропия Крылова-Колмогорова в методе молекулярной динамики является свойством системы и не зависит от параметров вычислительной схемы.
Рассматриваемая система - перемешивающая. В подобных системах возникает так называемый детерминированный хаос. Для этих систем характерны необратимость и непредсказуемость. Из свойства перемешивания вытекает эргодичность движения, т.е. средние по времени совпадают со средними по пространству. Такое поведение системы атомов обусловлено следующими причинами: нелинейностью системы дифференциальных уравнений (1), локальной неустойчивостью движения по Ляпунову (вследствие чего близкие фазовые траектории экспоненциально расходятся и сложным образом запутываются в ограниченной области фазового пространства), а также принципиальной невозможностью точной постановки начальных условий (2) для задачи Коши (эти начальные условия по существу задаются лишь указанием некоторого диапазона значений начальных импульсов и начальных координат, что может быть обусловлено различными причинами, в частности, соотношением неопределенности Гейзенберга). Непредсказуемость и необратимость в подобных задачах естественным образом вытекают из новой трактовки задачи (1), (2), которая уже по существу не является задачей Коши в классическом смысле.
Таким образом, задача отыскания траекторий отдельных атомов (1), (2) на достаточно больших временах нуждается в серьезной корректировке. Видно, что точное вычисление траекторий атомов по методу молекулярной динамики принципиально невозможно, даже если бы использовалась воображаемая «идеальная» разностная схема с шагом по времени, который стремится к нулю. С другой стороны, данное фундаментальное свойство оправдывает применение метода молекулярной динамики для определения макроскопических параметров системы. Для расчета наблюдаемых макроскопических величин путем статистического усреднения в методе молекулярной динамики необходимо найти решение за времена, много большие величины 1/К.
В данной работе моделировалось создание аморфного металла меди Си путем сверхбыстрого охлаждения расплава. Расчетный блок содержал 13500 атомов и представлял собой участок тонкой пленки. Взаимодействия атомов описывались с помощью парного потенциала Морза
и(г) = е(в-2а<г-^) - 2в-а(г-^у)
с известными параметрами е, а, а для меди [8]. По двум осям задавались периодические граничные
условия, по третьей - свободные. Первоначальная температура составляла 4500К. Такое большое значение температуры выбрано, чтобы минимизировать время получения расплава. Шаг по времени был равен 0,01 пс = 10-14 с.
Как показывает опыт расчетов, выполненных разными авторами, первоначально достаточно расставить атомы в виде правильной решетки Си (ГЦК кристалл), а начальные импульсы выбрать одинаковыми по модулю (значения импульсов должны соответствовать выбранной температуре) со случайным распределением их направлений, выбранных, конечно, так, чтобы полный импульс расчетного блока был равен нулю. Вследствие ляпуновской неустойчивости траекторий отдельных атомов в расплаве первоначально формируется случайным образом та или иная атомная структура, содержащая отдельные атомы и скопления атомов - кластеры, а также поры. Структурные единицы расплава хаотически перемещаются внутри расчетного блока. Время плавления - 10 пс. Характерное время расцепления корреляций 1/К составило величину порядка пикосекунды.
Численные расчеты показывают, что структурный хаос в жидком состоянии в значительной степени наследуется в твердом аморфном состоянии. При охлаждении расплава кинетическая энергия атомов быстро уменьшается и становится малой по сравнению с характерной потенциальной энергией их взаимодействия. При этом резко возрастает вязкость расплава, которую можно оценить на основе проведенных расчетов, а также экспериментальных данных. Наконец, вязкость расплава становится настолько большой, что структурные единицы расплава как бы «останавливаются» и у них сохраняются только внутренние степени свободы. При этом расплавленный металл затвердевает в аморфном состоянии, наследуя ту структуру, которой он обладал в жидком состоянии. Процесс кристаллизации становится невозможным, так как за очень короткое время атомы не успевают перемещаться на расстояния, которые позволили бы им образовать кристаллическую решетку. Такое метаста-бильное состояние не является состоянием с минимальной энергией, причем энергия образовавшегося аморфного металла определяется возникшей атомной структурой. Эти атомные структуры в значительной степени случайны, что обусловлено предшествующими им стохастическими процессами.
Далее аморфный металл охлаждался до низких температур, чтобы, в частности, исключить тепловые колебания при анализе атомной структуры. Численные расчеты демонстрируют возможность получения разнообразных атомных структур в зависимости от начальных условий и скорости охлаждения.
Структура расчетных блоков анализировалась с помощью специального визуализатора. Проводил-
ся анализ содержания фаз ГЦК, ГПУ и многогранников Франка-Каспера. Для каждого атома из расчетного блока оценивалось расположение ближайших соседей, после чего проводилось сравнение с эталонными образцами. Из фигур Франка-Каспера
больше всего обнаружено икосаэдров (16-вершинных фигур Франка-Каспера найдено не было).
В таблице приведены процентные содержания элементарных ячеек ГЦК, ГПУ и Франка-Каспера для трех различных скоростей охлаждения.
Процентные содержания элементарных ячеек ГЦК, ГПУ и Франка-Каспера для трех различных скоростей охлаждения
Скорость охлаждения Содержание элементарных ячеек Доля ГЦК среди элементарных ячеек Доля ГПУ среди элементарных ячеек Доля фигур Франка-Каспера среди элементарных ячеек
1013 К/с 29-79 2-80 0-59 7-94
1014 К/с 34-40 45-49 23-28 23-27
1015 К/с 3337 41-44 25-29 27-38
Процентное содержание - отношение числа атомов, являющихся «центрами» ячеек, к общему числу атомов в системе.
Важно отметить, что указанный в таблице разброс значений не является только следствием погрешностей расчетов, а скорее обусловлен случайным характером образования той или иной атомной структуры. Особенно большой разброс значений имеет место для скорости охлаждения 1013 К/с, когда образуется аморфно-нанокристаллическая структура. При такой скорости охлаждения вообще невозможно выделить доминирующую фазу, а случайный характер процесса образования атомных структур становится совершенно очевидным. Упорядоченная структура при этом занимает от 80 до 95% всего объема расчетного блока. Скорость охлаждения 1013 К/с близка к критической скорости охлаждения, условно разграничивающей области образования аморфной и кристаллической структур.
Необходимо отметить, что в соответствии с указанным выше случайным характером образующихся аморфных структур элементарные ячейки ГЦК, ГПУ и Франка-Каспера этих систем беспорядочно разбросаны по объему и не образуют сопряженных упорядоченных структур.
Как известно, выбор потенциала взаимодействия атомов играет важнейшую роль при применении метода молекулярной динамики (межчастичные силы можно в принципе рассчитывать путем приближенного решения уравнения Шредингера на
каждом временном шаге метода молекулярной динамики). Многими авторами указывается, что небольшие изменения потенциала могут значительно изменять результаты молекулярно-динамических расчетов. Если результаты при этом изменяются качественно, то рассматриваемую в режиме охлаждения совокупность атомов следует считать структурно неустойчивой (негрубой) системой (по Пон-трягину). Априори ясно, что структурную неустойчивость следует ожидать в окрестности критической скорости охлаждения. Полный теоретический анализ структурной устойчивости систем молекулярной динамики в настоящее время представляется трудно осуществимым. Однако качественная картина структурной устойчивости может быть получена путем численного анализа с вариацией эмпирических постоянных, входящих в потенциал межчас-тичного взаимодействия.
Таким образом, при исследовании атомных структур жидкостей и аморфных металлов необходимо использовать понятия и методы современной теории нелинейных динамических систем. На основании данных понятий и проведенных расчетов по методу молекулярной динамики приходим к выводу, что атомная структура аморфного металла принципиально непредсказуема как на уровне модели, так и в реальном эксперименте. Одновременно мы можем утверждать, что атомная структура аморфного металла принципиально невоспроизводима.
Авторы благодарят М.Д. Старостенкова за полезные обсуждения.
Библиографический список
1. Судзуки К., Фудзимори Х., Хасимото К. Аморфные металлы. - М., 1987.
2. Канн Р. У. Сплавы, быстрозакаленные из расплава // Физическое металловедение / под ред. Р. Кана. - М., 1987.
3. Краснов В.Ю., Полетаев Г.М., Старостенков М.Д. Исследование структуры аморфного никеля // Фундамен-
тальные проблемы современного материаловедения. -2006. - №4.
4. Квеглис Л.И. Структурообразование в аморфных и нанокристаллических пленках сплавов на основе переходных металлов: дис. ... д-ра физ.-мат. наук. - Красноярск, 2005.
5. Валуев А.А., Норман Г.Э., Подлипчук В.Ю. Метод молекулярной динамики: теория и приложения // Математическое моделирование: Физико-химические свойства вещества. - М., 1989.
6. Норман Г.Э., Стегайлов В.В. Стохастические свойства молекулярно-динамической ленард-джонсовской системы в равновесном и неравновесном состояниях // ЖЭТФ. - 2001. - Т. 119.
7. Заславский Г.М. Гамильтонов хаос и фрактальная динамика. - Ижевск, 2010.
8. Царегородцев А.И., Горлов Н.В., Демьянов Б.Ф., Старостенков М.Д. Атомная структура АФГ и ее влияние на состояние решетки вблизи дислокации в упорядоченных сплавах со сверхструктурой Ь12 // ФММ. - 1984. -Т. 58, №2.
