Содержание к диссертации
Введение
1 Обзор литературы 8
1.1 Метод молекулярной динамики 8
1.2 Равновесие фаз 16
1.3 Граница устойчивости фазы 17
1.4 Кинетика фазовых переходов первого рода 18
1.5 Крупномасштабное атомистическое моделирование фазовых превращений при внешних воздействиях 19
2 Гомогенный распад кристалла 22
2.1 Термодинамика метастабильных состояний 22
2.2 Статистическое описание 23
2.3 Метод расчета частоты гомогенной нуклеации 26
2.4 Результаты 33
2.5 Результаты главы 43
3 Стохастические свойства метода МД 45
3.1 Классификация приближений при решении МД задачи 46
3.2 Неустойчивость траекторий: времена вычислительной и динамической памяти 47
3.3 К-энтропия и время динамической памяти 53
3.4 Зависимость Kt^ от At и {АЕ2} 55
3.5 Выбор точности численного интегрирования 58
3.6 Физический смысл и роль времени динамической памяти 60
3.7 Время динамической памяти и характер предсказуемости времени жизни метастабильного состояния 62
3.8 Результаты главы 64
4 Распад и плавление перегретой кристаллической меди 69
4.1 Гомогенная нуклеация расплава в объеме 69
4.2 Температурная зависимость модулей упругости и условия устойчивости 75
4.3 Плавление с открытой поверхности 76
4.4 Результаты главы 83
5 Кавитации в жидком РЬ при отрицательных давлениях 85
5.1 Модель и метод расчета 86
5.2 Граница устойчивости метастабильной жидкой фазы 87
5.3 Частота кавитации 90
5.4 Обсуждение результатов 91
5.5 Результаты главы 95
Заключение 96
- Равновесие фаз
- Неустойчивость траекторий: времена вычислительной и динамической памяти
- Температурная зависимость модулей упругости и условия устойчивости
- Частота кавитации
Введение к работе
Диссертация посвящена разработке метода теоретического исследования динамики и кинетики фазовых превращений в метастабильных кристаллах и жидкостях на атомистическом уровне. Развитый подход примсттетт к неравновесным релаксационным процессам плавления перегретого кристалла и кавитации в растянутой жидкости. Исследованы стохастические свойства метода молекулярной динамики (МД), существенные для рассмотренных задач.
Актуальность работы. Развитие экспериментальной техники, освоение нано-, пико- и фемтосекундттых диапазоттов открывает возможности получения метастабильных состояний конденсированных веществ со значительными степенями перегрева и/или растяжения по сравнению с равновесным состоянием. Упомянем сильные ударные волны, нагрев фемтосскундным лазером, наносекундный электровзрыв проводников и др. Такое развитие науки ставит перед теорией новые задачи, например, о перегреве твердого тела с открытой поверхностью, который ранее в статистической физике считался, вообще говоря, невозможным. Требует более внимательного рассмотрения и классическая теория нуклеации (КТН), надежность которой может оказаться недостаточной, например, при импульсном растяжении жидкостей. Возникающие новые прикладные задачи также стимулируют развитие теории устойчивости метастабильных состояний и их распада.
Мощным инструментом развития теории конденсированного состояния является классический метод МД. Метод основан на решении классических уравнений движения многочастичной системы. Используя адекватные потенциалы межчастичного взаимодействия, можно исследовать на атомистическом уровтте широкий класс физических явлений в жидкостях, твердых телах, неидеалыгой плазме, биомолекулярных системах.
Важным потенциальным преимуществом метода МД является возможность исследования релаксационні.тх процессов в плотных средах исходя непосредственно из расчета динамики мпогочастичных систем без дополнительных предположений, присущих кинетической теории. Однако применение метода МД к изучению релаксационных процессов сравнительно менее развито по сравнению с изучением равновесных состояний. Не были даже сформулированы требования к МД расчетам релаксационных процессов.
Другим недостаточно изученным вопросом теории метода МД является соотношение динамических и стохастических свойств. Известно, что широкому классу динамических систем присущи хаотические свойства, в частности, экспоненциальная тгеустойчивость фазовых траекторий. По этой причине горизоттт предсказуемости во времени динамической эволюции мттогочастич-
j РОС НАЦИСЧМЛЬНАЯ і
3 | БИБЛИОТЕКА і
( С. Петербург /л \
*-оэ ти^г;
ной системы существенно ограничен. В связи с этим становится актуальным исследование характера предсказательных возможностей метода МД при исследовании релаксационных процессов и влияния на результаты конечной точности численных методов, использующихся при проведении расчетов.
Цель работы состоит в: 1) разработке метода МД исследования процессов распада метастабильных состояний с учетом результатов анализа стохастических свойств метода МД, 2) апробации разработанных подходов для расчета кинетики нуклеации п метастабильных кристаллах и жидкостях как для модельных систем, так и для металлов, 3) сравнении полученных результатов с имеющимися теоретическими подходами на основе КТН для описания кинетики фазовых превращений и 4) изучении фазовой диаграммы в области метастабильной кристаллической и жидкой фазы.
Научная новизна работы. Разработан молекулярно-динамический метод расчета частоты гомогенной нуклеации в метастабильной конденсированной фазе, основанный па формировании ансамбля независимых начальных условий и последующем усреднении времени жизни метастабильной фазы по ансамблю МД траекторий.
Проанализированы стохастические свойства МД. Показано, что метод МД является методом, который- і) сохраняет ньютоновскую динамику на времена меньше, чем время динамической памяти, и и) производит статистическое усреднение по начальным условиям вдоль МД траектории. Показано, что время динамической памяти растет лить логарифмически при увеличении точности расчета. Получены универсальные зависимости, связывающие энтропию Крылова-Колмогорова (усредненный максимальный показатель Ляпунова), время динамической памяти и уровень шумов в системе. Введено понятие времени вычислительной памяти — предельного горизонта предсказуемости решения конечно-разностной схемы.
С использованием разработанного подхода и реалистических многочастичных потенциалов межатомного взаимодействия проведен расчет частоты гомогенной нуклеации расплава и кавитации соответственно в перегретой кристаллической меди и растянутом расплаве свинца. Показано, что оценки по КТН могут быть приведены в соответствие с результатами МД расчетов.
Результаты расчетов температурной зависимости модулей упругости кристалла свидетельствуют, что граница устойчивости кристаллической меди определяется потерей устойчивости по отношению к сдвигу при повышении температуры. Показано соответствие механического и кинетического критериев устойчивости.
Оценены границы устойчивости перегретой кристаллической меди при положительных и расплава свинца при отрицательных давлениях.
Практическая ценность работы. Развитый подход может быть использован для определения частоты гомогенной и гетерогенной кавитации в жидкостях при отрицательных давлениях для использования в практически важных задачах гидродинамики. Уравнение состояния и границы устойчивости мстастабильных состояний могут быть включены в широкодиапазонные уравнения состояния вещества.
Положения, выносимые на защиту.
-
Метод расчета частоты гомогенной нуклеации. Разработанные процедуры усреднения.
-
Характер проявления стохастических свойств метода МД при расчете релаксации метастабильных состояний. Формула, связывающая время динамической памяти, усредненный максимальный показатель Ляпунова и уровень флуктуации полной энергии.
-
Частота гомогенной нуклеации в кристаллической меди, ее зависимость от степени перегрева. Частота гомогенной кавитации в расплаве свинца, её зависимость от степени растяжения.
-
Границы устойчивости перегретой кристаллической меди и растянутого расплава свинца.
Апробация работы. Основные результаты диссертационной работы были доложены на научно-технических конференциях МФТИ (Долгопрудный, 1999 - 2004), конференциях "Уравнения состояния вещества" (п. Эльбрус, 2002 и 2004), "Воздействие интенсивных потоков энергии на вещество" (п Эльбрус, 2003 и 2005), "Nucleation theory and applications" (Дубна, 2003, 2004, 2005), "Foundations of molecular modeling and simulation" (США, Keystone, 2003), "Фа-зовые превращения при высоких давлениях" (Черноголовка, 2004), Computational Physics (Германия, Aachen, 2001' США, San Diego, 2002; Италия, Genova, 2004), Computer Science (Нидерланды, Amsterdam, 2002), 15-th Symposium on thermophysical properties (США, Bolder, 2003), научно-координационных сессиях "Исследования неидеалытой плазмы" (Президиум РАН, Москва, 2001 и 2004), симпозиумах "Проблемы физики ультракоротких процессов в сильнонеравновесных средах" (Новый Афон, 2003, 2004, 2005).
Публикации. По материалам диссертации опубликовано 12 работ в реферируемых научных изданиях [1-12], 4 работы в сборниках [13-16] (список в конце автореферата) и тезисы российских и международных конференций.
Структура и объем диссертации. Диссертация состоит из введения, пяти глав и заключения, изложена на 111 страницах, включая 37 рисунков и 117 наименований цитируемой литературы.
Равновесие фаз
Условия равновесия при сосуществовании двух фаз состоит в равенстве соответствующих химических потенциалов обеих фаз при одинаковых давлении и температуре [31,36]. Метод молекулярной динамики позволяет рассчитать химический потенциал фазы путем термодинамического интегрирования (так называемого А-интегрирования [21]) с использованием начального состояния с известным абсолютным значением свободной энергии Гиббса. Например, кристалл Эйнштейновских осцилляторов и идеальный газ для кристаллического твердого тела и жидкости соответственно. Таким образом были рассчитаны параметры фазового равновесия кристалл - жидкость для систем Лепнард-Джонсовских атомов [37], металлов с использованием классических [38] и ab initio [39] потенциалов взаимодействия. Существуют методы термодинамического интегрирования на основе МД, позволяющие в рамках одного МД расчета на основе одной известной пары параметров (например, Тт, Рт) получать параметры кривой фазового равновесия, определяемого уравнением Клайперона-Клаузиуса [40]. С использованием МД возможно непосредственное моделирование сосуществования двух фаз в расчетной ячейке. Данный подход является менее вычислительно ресурсоемким и позволяет определять параметры фазового равновесия с хорошей точностью (см., например, [41,42]). В работе [43] таким образом исследовалось равновесие жидкости и пара в Леннард-Джонсовской системе с плоской поверхностью раздела. 1.3. Граница устойчивости фазы Для устойчивости однородной системы относительно малых (непрерывных) изменений параметров состояния должна быть положительной квадратичная форма д2и (где и — удельная внутренняя энергия), записанная через вариации энтропии и объема [44]: При усложнении системы появляются дополнительные термодинамические переменные и условия, аналогичные (1Л8) и (1.19). Для твердого тела в приращение внутренней энергии войдут члены вида cr rfe ., описывающие работу деформации, где 0 и 6 — составляющие тензора напряжений и тензора деформаций. Их связь дается обобщенным законом Гука через тензор модулей упругости твердого тела: 0"ifc = Ciklmelm- (1.20) Для кубической кристаллической решетки при гидростатической нагрузке Р условия устойчивости могут быть записаны в виде (используя сокращенную запись по Фойгту) [40,45-48]: Связь термодинамического критерия устойчивости с критерием Лин-демана рассматривалась в работе [49]. Одним из основных методов описания кинетики фазовых переходов первого рода является классическая теория нуклеации (КТН). Разработаны подходы описания кавитации [50,51], гомогенного [52-55] и гетерогенного плавления [56].
Исследование процесса нуклеации на основе метода молекулярной динамики позволяет уточнить основные предположения КТН и рассчитать неизвестные из эксперимента параметры. В работе [57] проведено моделирование процесса спонтанной кавитации в растянутой жидкости. По уменьшению до нуля времени ожидания зародыша определяется кинетическая граница устойчивости метастабиль-ной фазы, при этом не проводится статистического анализа результатов, время жизни метастабильной фазы оценивается по одиночной МД траектории. Проводится попытка сравнения результатов расчетов с классической теорией нуклеации. В работе [58] разработана теория, позволяющая рассчитатв зависимость поверхностного натяжения на границе пар-жидкость от степени кривизны поверхности. Определена зависимость работы зародышеобразования от размера зародыша. Результаты проведенного МД моделирования подтверждают выводы теории. В работе [59] максимальные степени перегрева и переохлаждения, достижимые при различных скоростях нагрева (охлаждения) исследованы на основе классической теории нуклеации, метода МД и динамических экспериментов. Показано, что максимально (минимально) достижимая температура перегретого твердого тела (или переохлажденной жидкости) зависит от безразмерного параметра, характеризующего нуклеационный барьер (который определяется по поверхностной энергии кристалл-расплав, энтальпии и температуре плавления), и скорости нагрева (охлаждения). Результаты МД моделирования и динамических экспериментов согласуются с предложенным подходом систематизации результатов. Современные вычислительные возможности позволяют проводить крупномасштабное МД моделирование систем, содержащих миллионы частиц в расчетной ячейке. Это позволяет исследовать задачи, традиционно относящиеся к области механики сплошной среды: распространение ударных волн и волн разрежения, рост трещин. Преимущество расчетов такого рода состоит в том, что в модель не приходится закладывать приближенных механизмов кинетики фазовых переходов, которые проходят "естественным" образом. В работе [60] метод МД использован для исследования структуры стационарной ударной волны, распространяющейся в направлениях (100), (110) и (111) в Ленпард-Джоисовском идеальном г.ц.к. кристалле. В отличие от ударных волн в газах и жидкостях профиль ударной волны в кристалле характеризуется осцилляциями локальной температуры, плотности и давления, что имеет место даже для ударных волн большой амплитуды, сопровождающихся плавлением.
Обнаружено существенное различие структуры фронта ударной волны и условий начала фазовых превращений для различных кристаллографических ориентации направления распространения ударной волны. В работе [61] представлены результаты МД моделирования разрушения тонкой пленки, нагретой фемтосекундиым лазерным импульсом. В рамках модели нагрев считается мгновенным, поскольку за время импульса не успевает произойти сколько-нибудь заметного смещения вещества. Расчет показал, что разрушение пленки происходит из-за взаимодействия волн разгрузки и может рассматриваться как модель более сложного процесса откола тонкого поверхностного слоя мишени в случае, когда после прогрева слой остается твердым. В частности было обнаружено, что разрушение кристаллического порядка, вызванное деформацией растяжения и сильной анизотропией остаточных напряжений, приводит к разделению отрывающегося от мишени слоя на две части. При этом растяжение и формирование анизотропных напряжений происходит вследствие расширения нагретого кристалла. В работе [62] в рамках крупномасштабного (1.2 млн. частиц) МД моделирования был изучен механизм образования дефектов при пластической деформации за фронтом ударной волны в Лен нард-Джонсовском монокристалле. Было обнаружено, что в узком слое за фронтом ударной волны (несколько периодов решетки) в результате тепловых флуктуации образуются маленькие дислокационные кольца, был определен их критический размер и проведена оценка энергии активации процесса. В работе [63] приводится обзор результатов по МД моделированию динамики разрушения твердых тел. Рассмотрены распространение трещин в хрупких (brittle) и ковких (ductile) материалах, взаимодействие дислокаций при механическом упрочнении, до- и сверхзвуковое распространение трещин. В работе [64] было проведено исследование кинетики и микроскопических механизмов плаления и разрушения тонких пленок Ni и Аи под действием коротких (150 - 200 фс) лазерных импульсов. На основе метода МД и методов описания сплошной среды была разработана модель для описания поглощения излучения электронами проводимости и нагрева решетки в результате электрон-фононного взаимодействия. В работе было показано, что при интенсивностях излучения близких к порогу плавления, структурная эволюция пленок зависит от соотношения двух процессов: распространения межфазных границ кристалл-жидкость с внешней поверхности пленки и гомогенной нуклеации и роста зародышей расплава внутри кристалла. При этом относительный вклад механизмов гомогенного и гетерогенного плавления зависит от интенсивности излучения, длительности импульса и параметров, характеризующих электрон-фононную релаксацию.
Неустойчивость траекторий: времена вычислительной и динамической памяти
Возникновение необратимости, стохастичности и хаоса в динамических системах исследовалось во многих работах, см., например, [19,20,31, 73-108] и содержащиеся там ссылки. Большой интерес при этом представляют расходимость траекторий и характеризующая эту расходимость К-эптропия (энтропия Крылова-Колмогорова, усредненный максимальный показатель Ляпунова). Величина К является также скоростью роста энтропии [73,80,81], т.е. /Г-1 оказывается важным временем релаксации. Было введено время предсказуемости поведения трг [75-77]. Это время характеризует интервал, на который можно предсказывать будущее поведение динамической системы, исходя из начальных условий и детерминистических динамических уравнений, определяющих эволюцию системы. В качестве причин конечности rw отмечались [75-78] измерительные шумы, флуктуационные силы и неточность знания дифференциальных уравнений, описывающих динамическую систему. Горизонтом предсказуемости 7 было названо значение, к которому стремится трг в пределе пренебрежимо малых измерительных шумов и шумов незнания. Как трг, так и Th пропорциональны А+1, где А+ — наибольший положительный показатель Ляпунова, а коэффициент пропорциональности логарифмически зависит от уровня шумов. Полагалось также, что время корреляции системы V 0-5А;1 [75-77]. Заславский [73] пользовался терминами К-эптропия и время расцепления корреляций т. Значение К соответствует величине А+, усредненной по фазовому пространству, а г - времени 7. При этом, однако, полагалось Частным случаем динамических систем являются классические системы многих частиц. Система уравнений (1.4) или (1.5) экспоненциально неустойчива для системы, состоящей больше чем из двух частиц, см., например, [19,20,73,80,81,92-94]. Параметром, определяющим степень неустойчивости, т.е. скорость разбегания первоначально близких фазовых траекторий, является К-энтропия. Значения К-энтропии рассчитывались методом МД для систем нейтральных [20,81-85,88-90] и заряженных частиц: двух- [96-98] и однокомпопептпой [99-101] плазмы. В [102] введено понятие о времени tfn динамической памяти, характеризующем время забывания начальных условий. Значение t6m определяется точностью схемы численного интегрирования [20,88,89,96-102]. Величины tfn и rh [75-77] имеют близкую природу, их различие обусловлено различием шумов, которым они соответствуют.
В данной главе рассматривается 3-х мерная системы частиц, взаимодействующих по потенциалу Леннарда-Джонса. Использовались периодические граничные условия. Число частиц в основной ячейке N менялось от 16 до 216. В рассчитывались и tfa, найдена зависимость величины от флуктуации полной энергии системы Е. Эта зависимость сопоставляется с результатами для плазмы. Обсуждается физический смысл К-энтропии. Физический смысл максимального динамического времени памяти в реальных системах связывается с малыми, но конечными шумами квантовой неопределенности, существующими в любой классической системе. Этот результат имеет отношение к гипотезе Ландау [31,74] о квантовой природе необратимости. Анализируется влияние конечной точности машинного представления чисел на возможности предсказуемости численного решения уравнений движения. Решение указанной системы дифференциальных уравнений может быть найдено лишь численно (начиная с N 3) на основе конечно-разностной аппроксимации (см., например, (1.14)). Заданный тип схемы и шаг интегрирования однозначно определяют функцию { (Atk), - (Atk) = О,1,...} (п — обозначает "численная"). В случае парных потенциалов сила действующая на г-ю частицу является суммой вкладов всех ее соседей jt=i где Gi — упорядоченная последовательность индексов. Структура Gi зависит от алгоритма сортировки частиц при расчете сил. При суммировании па компьютере с конечной точностью представления чисел (фиксированным числом десятичных знаков) различные поряд-ки суммирования в (3.1) дают различный результат F{. Обозначим траекторию полученную на заданном вычислительном устройстве {т (Аік), { (Atk), к = 0,1,...} (с — обозначает "расчетная"). Б силу Ляпуновской неустойчивости две МД траектории, рассчитанных из близких, но не совпадающих начальных условий, экспоненциально разбегаются со временем. Пусть {ri(t), Vi(t)) обозначает 1-ю траекторию и ( (), (0) обозначает 2-ю траекторию, рассчитанную из слабо возмущенной начальной конфибудут возрастать со временем экспоненциально: где ti некоторое переходное время (порядка обратной частоты межчастичных соударений), А, В и t[ зависят от степени возмущения начальной конфигурации ((Дг2(0)}, Дг;2(0))), К — значение К-энтропии.
Отметим, что значения К-энтропии слабо зависят от N [90], начиная с N 10, рис. 3.1. Через некоторый промежуток времени разбегание траекторий становится порядка характерного размера систе- — тепловая скорость. Характер разбегания при этом изменяется: {Ar2(t)\ переходит в дифузионный режим, a (Av2(t)) выходит на постоянное значение: где D — коэффициент диффузии. Время tm представляет собой время памяти, т.е. время, в течение которого МД система помнит тот факт, что начальные конфигурации на обеих траекториях были близки в момент времени t — 0. Так же как А и В, время памяти tm зависит от степени начального возмущения. Время вычислительной памяти. Ошибки численного округления являются слабым неконтролируемым источником возмущения при расчете МД траектории. Для исследования их влияния на результат расчета была использована процедура искусственной перестановки индексов при суммировании в (3.1). Путем случайных перестановок порядка суммирования в (3.1) реализуются различные варианты ошибок округления, и, следовательно, различные значения F;. Пусть две траектории {г01 (Atk), Vе1 (Atk), к — 0,1,...} и {гс2(Ді ),іЯ2(Д(А;), к = 0,1,...} рассчитываются из одной и той же начальной конфигурации, но порядок суммирования при вычислении сумм вида (3.1) разный. Траектории разбегаются экспоненциально, т.к. ошибки округления действуют как источники начального возмущения. Будем называть время памяти в этом случае временем вычислительной памяти tcm. Это время зависит от точности машинного представления действительных чисел, т.е., например, различным типам данных в языке программирования C++: float, double, long double (см. рис. 3.2). Величина tcm определяет время потери корреляции рассчитанной траектории { (Atk), v{Atk)} и точного решения конечно-разностного приближения {fn(Atk),vn(Atk)} для одной и той же начальной конфигурации.
Температурная зависимость модулей упругости и условия устойчивости
Был проведен расчет температурной зависимости модулей упругости при приближении к границе устойчивости кристалла (т — 0). Для расчет модулей упругости использовались флуктуационныс формулы с учетом температурного вклада и Борковских слагаемых согласно [110,111]. При этом усреднение значений флуктуационных и Борновских вкладов проводилось только по участку МД траектории до начала гомогенной нуклеации. Поэтому полученные значения модулей упругости можно отнести к метастабильному кристаллу как таковому и таким образом проверить выполнение условий устойчивости (1.21). Расчет показал, что модуль сдвига G = (Сц — С\ч — 2Р)/2 — 0 при т —» 0 (см. рис. 4.3). Таким образом, граница устойчивости кристалла определятся потерей устойчивости по отношению к сдвигу. Отметим, что при этом модуль всестороннего сжатия К = (Сц+2(7і2+-Р)/3 существенно отличен от нуля. Одновременное обращение в ноль модуля сдвига кристалла и среднего значения времени жизни говорит об эквивалентности термодинамического и кинетического критериев устойчивости перегретого кристалла, что а priori не является очевидным (см., например, [53]). На рис. 4.4 показаны зависимости величин К (Сц + 2С\г + -Р)/3, G = (Сц — С\2 — 2Р)/2 и G = С 44 — Р от температуры вдоль трех изохор, показанных в свою очередь на .РТ-диаграмме на рис. 4.5. Видно, что во всех трех случаях граница устойчивости определяется обращением в ноль G . На РТ-диаграмме показана полученная таким образом граница устойчивости кристаллической меди. Заметим, что в этом случае не имеет места геометрическое положение спииодали, как огибающей семейства изохор, справедливое для спииодали метастабильных жидкой и газовой фазы (см., например, [44]). В связи с тем, что кристаллическая решетка без внутренних дефектов и межзеренных границ устойчива по отношению к гомогенной нуклеации в широком диапазоне перегревов, определяющую роль при плавлении в экспериментальных условиях играет гетерогенное плавление. Вообще говоря, гетерогенное плавление может проходить практически безбарьерно, так как расплав обычно полностью смачивает поверхность твердой фазы и образование жидкости на неоднородностях начинается при любом перегреве ДГ 0. В импульсных процессах, когда интенсивности энерговклада в твердую фазу велики, нагрев кристаллической решетки может происходить на- столько быстро, что процессы гетерогенного плавления не успеют охватить значительную часть образца.
В связи с этим представляет интерес исследовать скорость распространения фронта плавления с неоднородностей, что делается в данной работе на примере расчета температурной зависимости скорости фронта плавления с открытой поверхности. Метод расчета основывается на работах [38,112]. Расчет выполняется с реалистическим многочастичным потенциалом для меди [28]. Для моделирования кристалла с открытой поверхностью использовалась расчетная ячейка, вытянутая по направлению z с периодическими граничными условиями в направлениях хну. В первую очередь определяются параметры нулевой изобары. Для этого при N — 2048 рассчитывается семейство изотерм P(V(a),T), Т — 300, 500, ...2300 К. При расчете каждой точки с использованием термостата система выводится на равновесие при заданной температуре в течении 5000 шагов, значение давления усреднялось по следующим 5000 шагам интегрирования. По семейству изотерм определяется зависимость постоянной решетки от температуры а а(Т) при Р = 0. Полученная зависимость позволяет определить коэффициент теплового расширения и модуль объемного сжатия кристалла при Т — 300 К: а = 17.8 К-1 и К — 136 ГПа (в хорошем согласии с экспериментальными значениями для твердой меди 16.8 К"1 и 142 ГПа). Для моделирования плавления с поверхности при фиксированном перегреве кристалла AT = Т — Тт генерировалась г.ц.к. решетка N — 5000 частиц с постоянной а(Т), вытянутая в пропорции 5:1:1 (50 х 10 х 10 элементарных ячеек). На протяжении 5000 шагов интегрирования решетка приводилась в равновесие при заданной температуре Т в 3-х мерных периодических граничных условиях. При этом использовались ограничения на движение частиц (в виде ячейки Вигнера-Зейтца) для предотвращения неконтролируемого спонтанного образования дефектов или нуклеации расплава (при больших значениях перегрева). Полученная конфигурация используется для расчета в NVT ансамбле в 2-х мерных периодических граничных условиях. Типичный результат расчета приводится на рис. 4.6. Для плавления кристалла при постоянной температуре необходимо совершить работу, равную энтальпии плавления AhsL- В расчете AhsL определяется по увеличению полной энергии системы (за счет термостата). Полученное значение АЕШ = 0.12 эВ/атом= 12 кДж/моль находится в хорошем согласии с экспериментальным значением для меди 13 кДж/моль. В приведенном случае степень перегрева достаточно высока для гомогенного образования расплава. При меньших перегревах образование расплава идет только на поверхности кристалла. Для количественного анализа движения фронта плавления расчетная ячейка делилась на слои (см. рис. 4.7), в каждом из которых вычислялся статический структурный фактор где (fj) — координаты частиц, усредненные по некоторому временному интервалу для сглаживания влияния тепловых флуктуации, а к — вектор обратной решетки (в расчетах использовался к \\ [111]). Определенный таким образом статический структурный фактор позволяет различать обла-сти кристаллической и жидкой фазы в ячейке: 3(к) — 1 для идеальной решетки { — Ri, S(k) — 0 в полностью неупорядоченном состоянии (пороговым значением считается S = 0.3). Иллюстрация результатов расчета движения фронта плавления приводится на рис. 4.8. Граница раздела кристаллической и жидкой фаз с течением времени распространяется вглубь ячейки. В течении процесса скорость продвижения фронта практически на изменяется, хотя при перегревах порядка 12 - 15 % картина становится менее регулярной.
Скорость фронта плавления неодинакова с разных сторон ячейки и при различных начальных конфигурациях Го, поэтому требуется статистическое усреднение по нескольким независимым реализациям процесса при одинаковом значении перегрева. Результаты расчета зависимости скорости фронта плавления от тем-псратуры представлены на рис. 4.9. Максимальное значение перегрева 200 К (15%), использованное в расчетах, обусловлено началом гомогенного образования расплава и большими флуктуациями скорости фронта плавления. В пределе малых перегревов зависимость скорости фронта плавления от перегрева Vfront AT. Линейная аппроксимация расчетных данных позволяет оценить равновесную температуру плавления vfront{T — 0) — 1366 ± 10 К, что хорошо согласуется с температурой плавления меди 1356.6 К. Сосуществование кристалла с расплавом Как было отмечено, при расчете зависимости скорости фронта плавления от температуры постоянная температура поддерживается при помощи термостата (согласно (1.14)). В этом случае в конце расчета система полностью переходит в жидкую фазу. При расчете в адиабатических условиях при соответственных значениях первоначального перегрева и размера системы может происходить релаксация к двухфазному состоянию, когда граница раздела кристалла и расплава приходит в равновесие. Значение температуры сосуществования представляет собой непосредственный способ оценки температуры плавления кристалла. На рис. 4.10 представлены результаты расчета. Температура сосуществования фаз Т — 1323 К несколько хуже согласуется с температурой плавления меди, чем значение, полученное в предыдущем разделе. Причиной этого может быть тот факт, что при релаксации температуры из-за влияния периодических граничных условий возникают напряжения в направлениях х и у. В этом случае в процедуре молекулярно-динамического расчета необходим учет возможности динамического изменения размеров расчетной ячейки. 4.4. Результаты главы Проведено моделирование гомогенной нуклеации расплава в перегретой кристаллической меди с использованием реалистического многочастичного потенциала для меди.
Частота кавитации
Кавитация в рассматриваемой модели имеет случайный характер гомогенного зародышеобразования, так как в ПГУ отсутствуют поверхностные эффекты и пространственные неоднородности. Момент начала образования полости зависит от локальных флуктуации скоростей частиц и расстояний между ними, которые и образуют зародыш новой фазы (полость). По визуальным оценкам в исследованном диапазоне значений плотности и температуры размер критического зародыша соответствует объему не более чем 10 — 100 атомов, поэтому ПГУ не влияют на начальную стадию кавитационного процесса. Взрывной рост полости приводит к быстрой разгрузке растянутой жидкости, при этом влияние ПГУ становится существенным и требует отдельного изучения. На конкретной МД траектории время жизни однородной метаста-билыюй жидкой фазы зависит от начальной конфигурации и начального распределения скоростей частиц, а при идентичных начальных условиях от шага интегрирования. Статистическое усреднение для заданного термодинамического состояния [р,Т,Р) проводится по ансамблю М независимых начальных микроконфигураций, для каждой из которых определяется соответствующее время жизни. На рис. 5.3 приведены примеры распределений, получающиеся в расчетах. Видно, что модель Пуассоновского случайного процесса достаточно хорошо описывает рассматриваемый процесс кавитации. Получение распределений с выраженной экспоненциальной формой требует набора большой статистики (М 100 — 200). Для показательного закона распределения среднеквадратичная погрешность определения f по М измерениям есть Of = f/VM. Скорость спонтанного фазового превращения принято характеризовать средним числом критических зародышей, образующихся в единице объема в единицу времени, т.е. частотой кавитации J, которая рассчитывается как J — (fV)_1. Результаты расчетов представлены на рис. 5.4. 5.4. Обсуждение результатов Представляет интерес сравнить результаты расчетов с предсказаниями КТН. Для сравнения используется зависимость скорости нуклеации от температуры, основанная на подходе Зельдовича (см. [50,51]): где a — поверхностное натяжение на линии испарения при температуре Г, W — работа образования критического зародыша, Р — давление пара в критическом зародыше. Приближение (5.1) выбрано из-за простоты пред-экспонснциального множителя, в который из специфических параметров, характеризующих метастабильпую жидкость, входит лишь а. При сравнении (5.1) с результатами МД расчетов считаем, что давление в системе, содержащей критический зародыш равно среднему давлению на метастабильном участке (0 t т), а давление пара в критическом пузырьке пренебрежимо мало Р С Р. Используются экспериментальные данные по поверхностному натяжению расплава свинца на линии испарения [117].
Поверхностное натяжение чрезвычайно сильно влияет на работу образования критического зародыша, а следовательно, и на зависимость частоты нуклеации от температуры и давления. Поэтому разброс экспериментальных данных учитывается в виде доверительного интервала значений amin а (ттах, эвристически построенного на основе неопределенности линейной аппроксимации экспериментальных данных в область высоких температур (см. вставку на рис. 5.4). Таким образом, на основе (5.1) для каждой температуры определяются области J(P; и), amin о ттах. Из рис. 5.4 следует, что результаты расчетов качественно верно согласуются с оценкой на основе КТН, при этом количественное согласие ухудшается с ростом температуры. В рамках подхода (5.1) различие может быть интерпретировано, как систематическое занижение работы образования критического зародыша W. Расхождение теории и результатов расчетов может быть уменьшено при учете зависимости поверхностного натяжения от степени кривизны поверхности (см., например, [50,51]). При этом результаты расчетов могут свидетельствовать о том, что поверхностное натяжение для пузырьков критического размера на 5-10% больше значения для случая плоской границы раздела фаз. Несмотря на приемлемое согласие результатов расчетов с рассмотренным приближением КТН, оценка разме- pa критического зародыша по КТН дает Nn = (р/т)(4тг/3)(2а/\Р\) w 1 атом (р/т = 2.82 х 1028 м"3, Т = 700К, а = 0.431 Н/м, Р = -3.89 ГПа), что, вообще говоря, свидетельствует о выходе за пределы применимости макроскопического подхода КТН в рассматриваемой области. Выполнено МД моделирование расплава свинца с использованием реалистического потенциала межатомного взаимодействия в рамках метода погруженного атома. Рассмотрены состояния при температурах меньше критической (Т О.ЬТср) и больших отрицательных давлениях. Рассчитана граница устойчивости жидкой фазы (спинодаль). На основе разработанной в работе методики проведен расчет частоты кавитации в растянутом метастабильном расплаве свинца. Проведено сравнение результатов МД расчетов с классической теорией нуклеации. Развиты теоретические основы исследования динамики и кинетики фазовых превращений в метастабильных кристаллах и жидкостях методом МД. Разработан метод расчета частоты гомогенной нуклеаци и/кавитации, основанный на формировании ансамбля независимых начальных конфигураций и последующем усреднении времени жизни мстастабиль-ной фазы по ансамблю МД траекторий ( 2.3). В качестве апробации метода получена зависимость частоты спонтанной гомогенной нуклеации от степени перегрева кристалла для случаев потенциалов межчастичного взаимодействия мягких сфер и Лен нард а-Джонса (рис. 2.5 и 2.6).
Проведено МД моделирование распада кристаллической решетки в условиях нагрева с постоянной скоростью. Проверена адекватность описания распада в рамках представления нуклеации как Пуассоновского случайного процесса (рис. 2.4 и 2.9). Проанализированы стохастические свойства метода МД и определены его возможности для исследования релаксационных процессов в плотных средах. Получены универсальные логарифмические зависимости, связывающие энтропию Крылова-Колмогорова, время динамической памяти tfn и уровень шумов в системе (формулы (3.14) и (3.15)). Метод МД: і) сохраняет ньютоновскую динамику на времена меньше, чем tfn, и ii) производит статистическое усреднение по начальным условиям вдоль МД траектории ( 3.6). Показано, что в силу принципиальной ограниченности величины tfn возможно лишь статистическое описание распада метастабильиого состояния ( 3.7). Рассмотрены гомогенная нуклеация и поверхностное плавление в перегретых кристаллах и кавитация в расплаве при отрицательных давлениях. Показано, что результаты расчетов кинетики гомогенной нуклеации в перегретой кристаллической меди и кавитации в растянутом расплаве свинца находятся в соответствии с оценками по КТН ( 4.1 и 5.4). Получена зависимость скорости фронта плавления от температуры (рис. 4.9). Найдена граница устойчивости перегретой кристаллической меди при О Р 15 ГПа (рис. 4.5): она определяется возникновением сдвиговой неустойчивости кристалла при повышении температуры. Показано соответствие механического и кинетического критериев устойчивости ( 4.2). Рассмотрен расплав свинца при температурах меньше критической (Т О.ЬТср) и больших отрицательных давлениях и оценена граница его устойчивости (рис. 5.1). 6.2. Достоверность результатов Проверялась достоверность полученных результатов. Одни и те же или сходные по смыслу величины рассчитывались разными способами: а) предельно достижимый перегрев Th кристалла в зависимости от скорости нагрева рассчитывался непосредственно и по формуле, связывающей величину Th с частотой спонтанной гомогенной нуклеации J(T), определенной в другом расчете ( 2.4, п.2); б) установлено совпадение границ устойчивости, найденных из механического и кинетического критериев ( 4.2).
