Теоретические основы исследования методом молекулярной динамики фазовых превращений в метастабильных кристаллах и жидкостях тема автореферата и диссертации по физике, 01.04.02 ВАК РФ

На правах рукописи УДК 536.75^538.9

Стегайлов Владимир Владимирович

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

Специальность 01.04.02 — Теоретическая физика

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Москва - 2005

Работа выполнена в Московском физико-техническом институте на кафедре физики высоких плотностей энергии (базовый институт ИТЭС ОИВТ РАН).

Научный руководитель: доктор физико-математических наук,

профессор Генри Эдгарович Норман

Официальные оппоненты: доктор технических наук,

профессор Давид Кириллович Белащенко

доктор физико-математических наук, профессор Александр Юрьевич Лоскутов

Ведущая организация: Институт теоретической физики

им. Л. Д. Ландау РАН

Защита состоится 22 декабря 2005 г. в 10 часов на заседании диссертационного совета К.212.156.05 при Московском физико-техническом институте по адресу: 141700 Московская обл., г. Долгопрудный, Институтский пер. 9.

С диссертацией можно ознакомиться в библиотеке МФТИ.

Автореферат разослан 17 ноября 2005 г.

Ученый секретарь диссертационного совета кандидат физико-математических наук, доцент

С. М. Коршунов

1 Общая характеристика работы

Диссертация посвящена разработке метода теоретического исследования динамики и кинетики фазовых превращений в метастабильных кристаллах и жидкостях на атомистическом уровне. Развитый подход применен к неравновесным релаксационным процессам плавления перегретого кристалла и кавитации в растянутой жидкости. Исследованы стохастические свойства метода молекулярной динамики (МД), существенные для рассмотренных задач.

Актуальность работы. Развитие экспериментальной техники, освоение нано-, ттико- и фемтосекундпых диапазонов открывает возможности получения метастабильных состояний конденсированных веществ со значительными степенями перегрева и/или растяжения по сравнению с равновесным состоянием. Упомянем сильные ударные волны, нагрев фемтосскундным лазером, наносекундный электровзрыв проводников и др. Такое развитие науки ставит перед теорией новые задачи, например, о перегреве твердого тела с открытой поверхностью, который ранее в статистической физике считался, вообще говоря, невозможным. Требует более внимательного рассмотрения и классическая теория нуклеации (КТН), надежность которой может оказаться недостаточной, например, при импульсном растяжении жидкостей. Возникающие новые прикладные задачи также стимулируют развитие теории устойчивости метастабильных состояний и их распада.

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

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

Другим недостаточно изученным вопросом теории метода МД является соотношение динамических и стохастических свойств. Известно, что широкому классу динамических систем присущи хаотические свойства, в частности, экспоненциальная неустойчивость фазовых траекторий. По этой причине горизоттт предсказуемости во времени динамической эволюции мттогочастич-

РОС НАЦИОНАЛЬНАЯ

БИБЛИОТЕКА

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

Цель работы состоит в: 1) разработке метода МД исследования процессов распада метастабильных состояний с учетом результатов анализа стохастических свойств метода МД, 2) апробации разработанных подходов для расчета кинетики нуклеации п метастабильных кристаллах и жидкостях как для модельных систем, так и для металлов, 3) сравнении полученных результатов с имеющимися теоретическими подходами на основе КТН для описания кинетики фазовых превращений и 4) изучении фазовой диаграммы в области метастабильной кристаллической и жидкой фазы.

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

Проанализированы стохастические свойства МД. Показано, что метод МД является методом, который- 1) сохраняет ньютоновскую динамику на времена меньше, чем время динамической памяти, и и) производит статистическое усреднение по начальным условиям вдоль МД траектории. Показано, что время динамической памяти растет лить логарифмически при увеличении точности расчета. Получены универсальные зависимости, связывающие энтропию Крылова-Колмогорова (усредненный максимальный показатель Ляпунова), время динамической памяти и уровень шумов в системе. Введено понятие времени вычислительной памяти — предельного горизонта предсказуемости решения конечно-разностной схемы.

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

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

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

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

Положения, выносимые на защиту.

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 наименований цитируемой литературы.

2 Содержание работы

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

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

Во второй главе описывается предлагаемый метод расчета частоты гомогенной пуклеации на примере нуклеации в перегретом кристалле и приводятся результаты для случая потенциалов мягких сфер и Леннарда-Джонса.

Метастабильное состояние можно считать термодинамически равновесным, если выполняются неравенства i, t* < f, где I, - характерное время релаксации по всем параметрам (температуре, давлению и др.), кроме возникновения новой фазы, t* — характерное время опыта, f - среднее время ожидания первого критического зародыша более устойчивой фазы в рассматриваемом объеме V. Распад перегретой метастабильной кристаллической решетки соответствует гомогенной пуклеации т.е. спонтанному флуктуацион-ному образованию зародышей жидкости и их росту. Процесс характеризуется частотой нуклеации J, которая равна среднему числу критических зародышей, образовавшихся в единицу времени в единице объема J = lj{fV).

В работе рассмотрена система N частит;, образующих ГЦК решетку при температуре Т в объеме V в 3-х мерных периодических граничных условиях. Для стабильных состояний кристалла при Т < Тт {Тгп — температура плавления) плотность вероятности в фазовом пространстве р„ (Г) соответствует каноническому NVT распределению и не меняется со временем (Г — точка в 6N-мерном фазовом пространстве). Распад метастабильного кристалла соответствует релаксации неравновесного распределения рП1 „ (Г) к равновесному ptiq (Г), соответствующему жидкой фазе. Эволюция рт п- (Г) —* ptiq (Г) может рассматриваться как движение отдельных фазовых точек Г = Г(£), Г(0) = Го. Каждой траектории Г(£) соответствует определенное значение времени жизни метастабильного кристалла г = т(Г0). Среднее время жизни f=/r(r)Ancr(r)rfT.

В работе предложена процедура аппроксимации рт, , (Г) при помощи со-

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

1. Идеальная рептетка N частиц генерируется для фиксированного значения постоянной решетки. Скорости частиц распределяются случайным образом. Система приводится к равновесию при Т < Тт.

2. Необходимая степень перегрева достигается при помощи включения термостата в схему расчета. Для предотвращения появления дефектов и сохранения решетки накладываются ограничения на движение частиц.

3. На фазовой траектории, полученной для заданной температуры в АТУТ ансамбле состояний, выбирается М независимых фазовых точек: Т30. ] = 1 ,..,М. Плотность вероятности ртсг( Г) = М-1 5(Г - Г^) где 5 — дельта функция, является аппроксимацией р„1СТ (Г). Точность аппроксимации возрастает с ростом М.

4. Для каждой выбранной фазовой точки Гц рассчитывается новая траектория, при этом отключается термостат, снимаются ограничения на движение частиц и находятся значения времени жизни т3.

5. В результате получается набор М значений времен жизни {ту}. Среднее время жизни по ансамблю есть г ~ + .. + тм)/М.

На рис. 1 (слева) показаны зависимости температуры и давления в МД системе от времени, характерные для распада метастабильного кристалла в расплав (результаты представлены в приведенных единицах параметров потенциала с и <т для энергии и длины и г = (гпсг2/с)1/'2 для времени, где тп — масса атомов). По этим зависимостям отчетливо видно, что на протяжении определенного времени кристалл находится в метастабильном перегретом состоянии, а затем спонтанно переходит в жидкую фазу. Времена жизни сильно различаются для разных начальных конфигураций Гд из одного и того же ансамбля для данного значения температуры. Полученный набор значений {т,} представляет собой результат М реализаций случайного процесса гомогенной нуклеации в перегретом кристалле. Соответствующую плотность вероятности можно оценить как ш^) ~ г??(£)/М = ^ -

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

Рис 1 • Слет. Зависимости мгновенных значений давления Р и температуры Т от времени па двух молекулярно-динамических траекториях (о и Д) Начальные конфигурации каждой из траекторий соответствуют одному ансамблю начальных конфигураций для температуры Т = 0 664 Скачок на приведенных зависимостях соответствует переходу метастабильного кристалла в жидкость. Времена жизни перегретого кристалла тта траекториях т, и т2 показаны стрелками. N — 6912 (потенциал мягких сфер, п 12). Справа: Распределения времен жизни т3 по ансамблю независимых начальных конфигураций (М = 30, N = 6912, потенциал Лепнарда-Джонса): а Т = 1 6097; б — Т = 1.5871. Метод построения распределений т (I) показан стрелками. Сплошные кривые построены по уравнению т(4) = (МП/т) ехр(—¿/г), где &1 — шаг разбиения (а <5г = 50; б — 100), г - среднее время жизни по ансамблю (а - г = 52; б - 133).

времени (/, 1 + сИ) не изменяется со временем и зависит только от параметров состояния.

По описанной методике были проведены расчеты зависимости частоты гомогенной нуклеации от температуры вдоль изохор для системы мягких сфер (рис. 2) и системы леннард-джонсовских частиц. Результаты расчетов частоты нуклеации 3 с различным числом частиц в молекулярно-динамической ячейке N при различных температурах ложатся на общую зависимость.

В расчетах с меньшим числом частиц достигаются большие степени перегрева. Анализ структурных превращений в расчетной МД ячейке показывает, что для N < (1 -г 3) х 103 распад проходит как одновременный коллапс всей решетки. При больших N различается момент возникновения зародыша жидкости и его рост. По визуальным оценкам критический зародыш жидкой фазы содержит несколько сотен частиц.

Образование жидкой фазы в перегретой кристаллической решетке связано с преодолением энергетического барьера XV, который определяется рабо-

J

с к

10"4

10'5

10-6

10-7

0.63 0.65 0.67 0.69 0.71 0.73 0.75

той образования зародыша критического размера. По данным МД расчетов может быть выполнена оценка значения числа Гиббса G = W/kgT для нук-леации в исследуемом диапазоне перегревов. Для системы мягких сфер в температурном диапазоне Т = 0.699 -г 0.735 значение G изменяется от 85 до 20. Для системы частиц, взамодействующих по потенциалу Леннарда-Джонся, оттенка среднего значения в интервале Т = 1.54 ~ 1.81 дает G ~ 45, что согласуется с оценкой G = 69 -г 22 на основе КТН по имеющимся данным для кристаллического аргона, взятым из работ Байдакова, Галашева, Скрипова (1980) и Скрипова, Файзуллина (2003).

Результаты для потенциала мягких сфер и потенциала Леннарда-Джонса показывают, что степени достижимых перегревов твердого тела по порядку величины составляют 20% —30%. Этот результат хорошо согласуется с имеющимися в литературе оценками (Lu, Li, 1998; Swift et al., 2003). Достижимая степень перегрева кристалла увеличивается с ростом давления.

Наряду со спонтанным распадом в работе изучался распад кристаллической решетки при нагреве. Использовалась схема расчета, которая также состояла из создания ансамбля начальных конфигураций и последующего усреднения по независимым МД траекториям. Температура кристалла увеличивалась с постоянной скоростью Т = dT/dt начиная с исходной температуры Т <Тт. При используемых в расчетах скоростях нагрева кристаллическая решетка находилась в состоянии термодинамического равновесия в каждый момент времени, что проверялось по соответствию распределения

Рис. 3' Слева: Зависимость мгновенной температуры и давления от времени на двух МД траекториях, рассчитанных при Т = 1 25 х 10-3 из различных начальных конфигураций. Показаны моменты начала распада Т'/1 и Т2Л (потенциал мягких сфер, п — 12). Справа■ Распределения показывают число МД траекторий из ансамбля М траекторий, на которых распад начинается в температурном интервале (Т. Т+АТ), АТ = 0.03 N = 4000: 1 — Т 6 25 х Ю-4 (М = 20), 2 — Т- 2 х Ю-2 (М = 22). Сплошные кривые рассчитаны как МАТ^(Т) (см. (1))

скоростей атомов Максвелловскому.

Выли получены распределения значений температуры Т'\ которая может быть достигнута до начала распада кристаллической решетки, т.е. до появления первого критического зародыша (рис. 3, справа). Проверено соответствие полученных распределений модели пуклеации как Пуассоновского случайного процесса. В условиях неизотермического процесса в данном случае плотность вероятности появления первого критического зародыша как функция температуры в однородной системе имеет вид (Скрипов, Коверда, 1984):

И(Г) = ехр ( - I Г J{T)dT). (1)

Т Г Jт0

Для проверки использовались результаты расчетов зависимости ./(Т). Результаты, полученные независимыми путями согласуются друг с другом.

В третьей главе проанализированы стохастические свойства метода МД.

С математической точки зрения метод МД заключается в решении системы обыкновенных дифференциальных уравнений при заданных начальных условиях (г°,у°) (точка в 6Д''-мерном фазовом пространстве), где =

(f^, ...,fдг) and v" = (i7J\ ... В практически важных случаях силы,

действующие на частицы, представляют собой гладкие функции координат, при этом справедлива теорема существования и единственности решения задачи Коши, т.е. существует точное решение (траектория) {^(t), ^(t)} для заданной начальной конфигурации (г°,гР).

Решение указанной системы дифференциальных уравнений может быть найдено лишь численно (начиная с N > 3) на основе консчтто-разностной аппроксимации. Заданный тип схемы и шаг интегрирования At однозначно определяют решение [r(t),v(t),t = к At, к = 0,1,...}. отличающееся от {^(i), из-за погрешностей конечно-разностной схемы.

Ляпуиовстя неустойчивость приводит к экспоненциальному разбеганию со временем двух МД траекторий, рассчитанных из близких, но не совпадающих начальных условий с одним и тем же шагом At, или рассчитанных из одних и тех же начальных условий, но разными шагами численного интегрирования At и At', что приводит к росту погрешностей численного решения уравнений движения. Пусть {rt(t). vt(t)} и обозначают 1-ю и 2-

ю траектории, Усреднённые отличия координат и скоростей первой и второй траектории в совпадающие моменты времени

{Агщ = jjf^im-кт <м*)> = ^¿(«w

г—1 г—1

возрастают со временехЧ экспоненциально (рис. 4, слсва)-

(Ar2(t)) = Aexp(Kt), {Av2{t)) = Bexp(Kt), ti<t< tm (3) где /,; некоторое переходное время, А, В и зависят от степени возмущения начальной конфигурации на нескольких первых шагах интегрирования, а величина К представляет собой усредненный по фазовому пространству максимальный показатель Ляпунова многочастичной системы (другие названия — энтропия Крылова-Колмогорова, К-энтропия). Значения К слабо зависят от N для N > 10 (Валуев, Норман, Подлипчук, 1990).

Через некоторый промежуток времени tm « K~l \п(В~^6квТ/т), когда разбегание достигает (Дг2(£)} ~ р-2/3 и (Av2(t)) ~ v2h, где р — N/V — плотность, vui = \/3кцТ/rri — тепловая скорость, характер разбегания изменяется: (Дг2(<)) переходит в дифузионный режим, а (Дv2(t)} выходит на постоянное значение, т.е. (Ar2(t)) = bD(t — tm) + (Ar2{tm)), (Дг>2{t)) = 2 (v2) = 2Vth, где D — коэффициент диффузии. Время tm представляет собой время памяти, т.е. время, в течение которого МД система помнит, что начальные конфигурации на обеих траекториях были одинаковы.

<AV2>

d«"

4Л <Ar?>

t

10-ю 10-20 10-30 1040 10-50 j ю-«"

1/ / / ■•' 2/ /

7/3

• / 4 /

• • /»/ *

/

<Ar2>

20 40 60 80 100

Рис. 4: Слева: Нормированные усредненные разбегания координат (Дг2(/)) и скоростей (Av2(t)) на двух траекториях, рассчитанных из тождественных начальных условий с шагами Д£=0.001 и Ai'=0.0001 (TV = 64, р = 0.5, Т = О 44, потенциал Леннарда-Джонса). Справа: Разбегания координат (Аг2(/)) на двух траекториях, рассчитанных из одних и тех же начальных условий при разном порядке суммирования F, Приводятся данные для различных вариантов машинного представления действительных чисел (N = 500, р = 0 45, Т — 2, потенциал Леннарда-Джонса с радиусом обрезания гс _ 21/б): 1 — 4 байта (PC), 2-8 байт (Cray SV1), S - 8 байт (PC), ^ - 16 байт (Cray SV1). Для случая наибольшей точности представления чисел (4) пунктиром показаны максимальное и минимальное значение разбегания траекторий отдельных частиц в даппый момент времени max, Arf(t) (1 =■ 1 ,..N).

Время динамической памяти tdm характеризует промежуток времени, за который теряется корреляция решения конечно-разностной аппроксимации {r(t),v(t)} и точного решения задачи Копти для той же начальной конфигурации {^(i),^^)}. Для определения времени динамической памяти проводятся расчёты tm при одном и том же значении At и различных значениях At': At/2, At/Ь, Ai/10 и т.д. (рис. 5). Предельное значение tm при (At/At') —> 0 есть время динамической памяти tfn для данной схемы и выбранного шага численного интегрирования At. В процессе численного интегрирования по прошествии времени t'jn система частиц полностью "забывает" свои начальные условия, и рассчитываемая МД траектория теряет всякую корреляцию с исходной ньютоновской траекторией.

Величина t'lm логарифмически растет при уменьшении шага численного интегрирования. Такой результат можно получить, предполагая, что А ~ (At)71, где п — определяется порядком точности схемы численного интегрирования. Действительно в момент времени t = &knT/m = 2(v2} =

* -

30

20

10

At'/At

Ktd„

♦ч

AE/E

о 1

0.2 0.3 0.4 0.5 10-7 10"6 10"5 10"4 10"3 Ю"2 10"1

Рис. 5. Слева: Зависимость tm от At'/At для различных значений At (N = 64, р - 0 5, Т - 0.44, потенциал Леннарда-Джонса): 1 — At - 0 01, 2 - 0.005, 3 — 0 001, 4 — 0.0005, 5 — 0.0001 (крайние левые точки на линиях 1 и 2 соответствуют At! - 0.00001). Справа• Зависимость безразмерного произведения Ktfn от относительной флуктуации полной энергии системы: кружки - Леннард-Джонсовская система (схема 2-го порядка точности), кресты - неидеальная плазма, ромбы - Леннард-Джонсовская система (схема 4-го порядка точности). Пунктир соответствует (4).

(Av2(tfn)) = Л cxp(Ktfn). Логарифмируя, получаем Ktfn = —nln( At) + const.

Вследствие приближённости численного интегрирования постоянство полной энергии системы Е имеет место лишь в среднем. Значение Е от шага к шагу флуктуирует около среднего значения, поэтому траектория, рассчитываемая в методе МД, лежит не на поверхности Е = const, как должно было бьт быть для решения уравнений Ньютона, а располагается в некотором слое толщины АЕ > 0 около поверхности Е = const (Валуев, Норман, Подлипчук, 1989). Значение АЕ определяется точностью и схемой численного интегрирования, причем (АЕ)2 ~ Atn. Следовательно,

Ktdm = -\n{(AE)2) + const. (4)

На рис. 5 (справа) представлены результаты для Леннард-Джонсовской системы, а также результаты для нсидеальной плазмы (Морозов, Норман, Валуев, 2001). Результаты расчетов хорошо согласуются с (4).

Время вычислительной памяти 1ст вводится в работе для характеристики промежутка времени, за который из-за ошибок округления теряется корреляция численного {r{t),v(L)} и точного {r"(£), v"(£)} решения конечно-разностной аппроксимации для одних и тех же начальных условий. Для

определения времени вычислительной памяти была использована процедура искусственной перестановки индексов: в случае парных потенциалов сила действующая на ?-ю частицу является суммой вкладов всех ее соседей

= YTk=\ Рчк'Зк 6 <Я> где С, — упорядоченная последовательность индексов. Структура 6'г зависит от алгоритма сортировки частиц при расчете сил. В результате конечной точности представления действительнт.тх чисел различит,ге порядки суммирования дают различный резулт.тат Рг. Различные варианты ошибок округления и, следовательно, различные значения реализуются путем перестановок порядка суммирования.

МД траектории, рассчитанные из одной и той же начальной конфигурации с одинаковым шагом численного интегрирования, но с разным порядком суммирования при расчете Ь\, разбегаются экспоненциально, т.к. ошибки округления действуют как источники возмущения (рис. 4, справа). Величина 1,ст зависит от машинного представления действительных чисел и для заданного представления определяет верхнюю границу возможности увеличения времени динамической памяти при повышении точности численттого интегрирования (уменьшении ДI и Д Е).

Статистический смысл метода МД. Из предыдущего следует, что время динамической памяти оказывается много меньше характерных времен МД расчета. В работе дается истолкование фактического смысла метода МД.

Сопоставим значение времени динамической памяти с характерными временами автокорреляционной функции скорости (АКС). Результаты для последней представлены на рис. 6 (слева). Область, когда значения АКС превышают Ю-1, соответствует временам, меньшим времени памяти Таким образом, корреляции в этой области являются динамическими корреляциями, следующими из уравнений Ньютона. Корреляции для хвоста АКС попадают в ту временную область, когда динамическая память о начальных условиях уже не сохраняется, т.е. эти корреляции имеют уже не динамическую, а стохастическую природу.

Конечность времени динамической памяти накладывает ограничения на горизонт предсказуемости во времени динамической эволюции многочастичной системы.

Таким образом, метод МД: 1) сохраняет ньютоновскую динамику на времена меньше, чем и п) производит статистическое усреднение по начальным условиям вдоль МД траектории.

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

Рис. 6 Слева: Автокорреляционная функция скорости (/V = 64, р — 0 5, Т - 0.44, потенциал Л еннарда^ Джонса).

Справа■ Распределения времен жизни перегретой кристаллической ГЦК решетки, полученные разными методами усреднения (Л' - 500, Т = 0 787, р - 1 28, потенциал мягких сфер, га — 8)' 1-е использованием ансамбля независимых начальных конфигураций, потученных при помощи расчета с удержанием атомов в сферических ячейках. 2-е использованием различных шагов численного интегрирования: Дг — 0 001, 0 0011, ... 0.0049, 0.0050. В обоих случаях число независимых траекторий М = 41. 3 — экспоненциальное распределение по уравнению т(т) — (АШ/т)ехр(-г/т), где Л мелкость шага распределения, т — среднее время жизни по ансамблю (г = 8 57).

степени перегрева. Только в непосредственной близости от границы устойчивости кристалла т ~ а для меньших степеней метастабилытости <С т. Таким образом, в большинстве случаев в методе МД не существует единственного времени жизни метастабильного состояния т(Го) для выбранных начальных условий Г0, следовавшего бы из классических детерменированных уравнений движения многочастичной системы.

Расчет показал, что при использовании различных шагов численного интегрирования разброс значений времени жизни может быть более одного порядка. Действительно МД траектории, рассчитываемые из Го с различными шагами численного интегрирования становятся независимыми за времена порядка Поэтому распределение значений времени жизни, полученное таким образом, должно соответствовать распределению, полученному для той же температуры путем усреднения по ансамблю независимых начальных условий (глава 2), и, в частности, соответствовать экспоненциальному Пуассоттовскому распределению. Расчет подтвердил справедливость этого предположения (рис. 6, справа).

Таким образом, метод МД дает возможно лишь статического описания релаксационных процессов на временах больших При этом практические возможности увеличения ограниченны: имеющиеся вычислительные возможности позволяют уменьшить АЕ только на 5 порядков даже при использовании рафинированных численных схем, что увеличило бы лишь вдвое, а кроме того и) конечная точность машинного представления действительных чисел приводит к ограничению < В силу принципиальной ограниченности величины возможно лишь статистическое описание распада метастабилыгого состояния. Вместе с тем задача чрезмерного увеличения не имеет значения для физики в силу приближенного характера классического описания динамики многоатомных систем.

Результаты рассмотренного примера характера предсказуемости при МД расчете распада метастабилыгого кристалла, вообще говоря, можно обобщить и на случай любых релаксационных процессов в плотных средах и выделить два метода усреднения. I) Усреднение по ансамблю независимых начальных условий. При этом процедура формирования такого ансамбля, соответствующего заданному неравновесному (метастабильному) состоянию представляет собой нетривиальную задачу. II) Усреднение по ансамблю траекторий, рассчитанных с различным шагом численного интегрирования, что может быть удобно в случае, если формирование ансамбля независимых начальных конфигураций затруднено.

В четвёртой главе предложенный метод расчета частоты гомогенной нуклеации используется для перегретой кристаллической меди. Модель соответствует г.ц.к. кристаллу N = 6912 частиц в 3-х мерных периодических граничных условиях. Используется многочастичный потенциал межатомного взаимодействия (М18Ып с! а!., 2001). Выполнены расчеты среднего времени жизни перегретой кристаллической решетки и частоты гомогенной нуклеации. Рис. 7 иллюстрирует процесс нуклеации на МД траектории.

Сравнение результатов МД расчета среднего времени жизни т = (./V)-1 и теоретических оценок на основе КТН представлено на рис. 8. Мотори-пым и Мушером (1984) были получено следующее выражение для скорости гомогенной нуклеации в твердом теле в случае больших перегревов:

л = шЛЩъ?), (6)

где п$ — концентрация атомов в твердой фазе, д — теплота плавления па единицу объема, ■>] — константа, характеризующая скорость распространения

!

!

Рис. 7: Спонтанный распад перегретого кристалла, температура Т = 1610 К (перегрев AT — 252 К, 18.6%). Показаны зависимости от времени мгновенного значения давления F и параметра порядка — статического структурного фактора |S(fc)|2. Появление и рост критического зародыша расплава проиллюстрированы на четырех ортогональных проекциях структуры частиц в расчетной ячейке в моменты времени, показанные стрелками (цвет частицы соответствует степени ее отклонения из соответствующего узла идеальной гц.к. решетки Дг,, а г„„ - расстояние между ближайшими соседями в идеальной решетке). Значение времени жизни перегретого кристалла на траектории 1 г — 165 пс; серым цветом показана другая траектория 2, начальная конфигурация которой слабо отличалась от начальной конфигурации 1- ой траектории. Флуктуацию в районе t = 70 — 80 пс можно трактовать как образование околокритического зародыша, который в случае 1 пе перешел через критический размер и исчез, а в случае 2 достиг критического размера

межфазной границы кристалл-расплав, а — энергия поверхностного натяжения межфазной границы кристалл-расплав. Энергия активации процесса гомогенного зародытнеобразования W — А/(к%(Т — Т„,)2), выраженной через а, существенным образом зависит от значения а, которое традиционно (из-за отсутствия дополнительных данных) принимается равным макроскопическому поверхностному натяжеиию ira границе кристалл-расплав при температуре Тт. Сделанная в указанных предположениях в работе Rethfeld'a с соавторами (2002) оценка температурной зависимости f (Т), показанная на рис. 8 (слева) пунктиром, качественно соответствует результатам МД расчетов, но

Рис. 8. Слева• Кружки — среднер время жизни перегретой кристаллической решетки т, V = 90 771 нм1. Пунктир — расчет по (5) и (б), сплошная линия - расчет по (5) и (7). Ромбы (штрих-пунктир) — результаты расчета модуля сдвига кристалла С. Справа ■ РТ-диаграмма Си. ЛИ — кривая равновесного плавления. Изохоры: 1 — а — 3 632А, 2 — а - 3 689А, 3 — а — 3 745А (где а — постоянная кристаллической решетки) CD — граница граница устойчивости кристаллической меди.

дает завышенное значение предельно допустимого перегрева (при f —» 0).

Однако, как было показано Моториным и Мушером (1984), на основе имеющихся результатов опытов (Коверда, Скрипов, 1973) по кристаллизации в переохлажденной меди при То ~ 0.8Тт можно определить температурную зависимость А(Т), пользуясь тем, что выражение (5) остается верным и для частоты гомогенной кристаллизации в расплаве. Тогда:

А(Т) = А(Т0) + (ЭА/дТ)\То(Т - То). (7)

Оценка т(Т), сделанная с учетом температурной зависимости А(Т) и показанная на рис. 8 (слева) сплошной линией, хорошо согласуется с результатами МД расчетов. Это говорит о согласованности описания кинетики гомогенной пуклеации в перегретой твердой меди методом МД (метод расчета и потенциал взаимодействия) и КТН с привлечением экспериментальных данных по кинетике обратного процесса — кристаллизации из расплава.

Был проведен расчет температурной зависимости модулей упругости при приближении к границе устойчивости кристалла (т —> 0). Использовались флуктуационные формулы с учетом температурного вклада и Борновских слагаемых (Ray, 1988; Wolf et al., 1992). Усреднение флуктуационных и Борновских вкладов проводилось по стационарному участку МД траектории до появления в системе критического зародыша. Расчет показал, что тетрагональный модуль сдвига G' —> 0 при f —> 0 (рис. 8, слева), т.е. граница устой-

чивости кристалла определятся потерей устойчивости по отношению к сдвигу. При этом модуль всестороннего сжатия К = -V{dP/8V)r существенно больше от нуля. Одновременное обращение в ноль модуля сдвига кристалла и среднего значения времени жизни говорит об эквивалентности механического и кинетического критериев устойчивости перегретого кристалла.

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

Для используемой модели меди рассчитана температурная зависимость скорости фронта плавления. Использован метод (Nguyen et al., 1992).

В пятой главе предложенный метод расчета частоты гомогенной нукле-ации используется для определения частоты гомогенной кавитации в растянутом расплаве свинца. Рассматривается система N = 13500 атомов в кубической расчетной ячейке в 3-х мерных периодических граничных условиях. Используется многочастичттый потенциал межатомного взаимодействия (Lim et al., 1992). Выполнены расчеты среднего времени жизни растянутой жидкой фазы и частоты гомогенной нуклеации. Путем использования различных исходных вариантов распределения скоростей формируется ансамбль независимых начальных условий, соответствующих жидкости при заданных р и Т. В исследованном диапазоне значений плотности и температуры размер критического зародыша соответствует объему не более, чем 10 — 100 атомов.

Для сравнения результаты расчетов с предсказаниями КТН используется зависимость скорости нуклеации от температуры, основанная на подходе Зельдовича (1942):

где а — поверхностное натяжение на линии испарения при температуре Т, Р' — давление пара в критическом зародыше. При сравнении (8) с результатами МД расчетов давление в системе, содержащей критический зародыш, приравнивалось среднему давлению па метастабильном участке (0 < £ < т), а давлением пара в критическом пузырьке пренебрегалось Р' Р. Значение поверхностного натяжения чрезвычайно сильно влияет па работу образования критического зародыша, поэтому разброс экспериментальных данных по поверхностному натяжению расплава свинца на линии испарения учитывается в виде доверительного интервала значений атт < а < <т„шх. на основе которого для каждой температуры по (8) определяются области 3(Р\о).

Из рис. 9 (справа) следует, что результаты расчетов качественно согласу-

(8)

Р,Ща

-0 5

-f

2000

4000 / j 6000 з 10» / '

/

/

/ 10»

Р,ГПау - -3

Т-1(РК /

/ -4 5

■5 /

3

-4

2

— p.i/i.M1

92 94 96 98

3-10»

102*

« i » »

;J,cmv4 \

4

t

i

\

т1 11

\ ■

\' •

\'\ л

f7

C7,MH/MV 3

Tjfcv

440 >

400-a

360-

32ol ,

Р,ГДа '

800 1200 1600 2000

-2 5

-15 -1

Риг. 9: Слева: PT-диаграмма Pb: 1 — экспериментальная кривая плавления и ее экстраполяция в область отрицательных давлений по уравнению Симона, 2 — кривая испарения, 3 — оценка положения спинодали жидкого свинца путем перенормировки спинодали Леннард-Джонсовской системы по параметрам критической точки CP (Тер — 5400К, Рср = 0 175ГПа, Терновой и др., 1996). Используется разрыв вертикальной оси и разные масштабы Показаны полученные в дапной работе точки спинодали (о) и состояния, в которых была рассчитана частота кавитации (•). На вставке приведен пример изотермы жидкого свинца (Т = 1000К). Сплошными линиями показаны экстраполяции МД результатов (♦) полиномами 2-5 степени Соответствующая неопределенность потоже-ния точки спинодали показана на РТ-диаграмме.

Справа: Зависимость частоты кавитации J от давления Р вдоль трех изотерм: 1 Т — 700К; 2 — 1000К; — 2000К. Результаты приводятся с погрешностью, соответствующей погрешности определения среднего времени жизни. Пунктиром отмечены границы областей па J — Р плоскости, которые соответствуют расчетам согласно (8) с учетом неопределенности значений поверхностного натяжения жидкого свинца а. На вставке показаны справочные значения <т: пупктиром изображена область значений, определенная по погрешности линейной экстраполяции экспериментальных точек; стрелки соответствуют неопределенности значений а при фиксированных значениях температуры.

ются с оценкой на основе КТН, при этом количественное согласие ухудшается с ростом температуры. Различие может быть интерпретировано, как занижение подходом (8) работы образования критического зародыша Оценка размера критического зародыша по КТН дает = (р/т)(4тт/3)(2сг/\Р\)3 « 1 атом, что также свидетельствует о выходе за пределы применимости макроскопического подхода КТН в рассматриваемой области.

Проведена оценка положения границы устойчивости (спинодали) в рамках рассматриваемой МД модели на основе расчета Р—р зависимостей вдоль изотерм и определения точки (дР/др)^ — 0 путем экстраполяции (см. вставку на рис. 9, слева).

3 Основные результаты и выводы работы

Развиты теоретические основы исследования динамики и кинетики фазовых превращений в метастабильных кристаллах и жидкостях методом МД.

• Разработан метод расчета частоты гомогенной пуклеации/кавитации, основанный на формировании ансамбля независимых начальных условий и последующем усреднении времени жизни метастабильной фазы по ансамблю МД траекторий. В качестве апробации метода получена зависимость частоты спонтанной гомогенной нуклеации от степени перегрева кристалла для случаев потенциалов межчастичного взаимодействия мягких сфер и Ленпарда-Джонса. Проведено МД моделирование распада кристаллической решетки в условиях нагрева с постоянной скоростью. Проверена адекватность описания распада в рамках представления нуклеации как Пуассоновского случайного процесса.

• Проанализированы стохастические свойства метода МД и определены его возможности для исследования релаксационных процессов в плотных средах. Получены универсальные логарифмические зависимости, связывающие энтропию Крылова-Колмогорова, время динамической памяти и уровень тумов в системе. Метод МД: 1) сохраняет ньютоновскую динамику на времена меньше, чем и п) производит статистическое усреднение по начальным условиям вдоль МД траектории. Показано, что в силу принципиальной ограниченности величины возможно лишь статистическое описание распада метастабильного состояния.

• Рассмотрены гомогенная нуклеация и поверхностное плавление в перегретых кристаллах и кавитация в расплаве при отрицательных давлениях. Показано, что результаты расчетов кинетики гомогенной нуклеации в перегретой кристаллической меди и кавитации в растянутом расплаве свинца находятся в соответствии с оценками по КТН. Получена зависимость скорости фронта плавления от степени перегрева.

• Найдена граница устойчивости перегретой кристаллической меди при О < Р < 15 ГПа: она определяется потерей устойчивости кристалла по отношению к сдвигу при повышении температуры. Показано соответствие механического и кинетического критериев устойчивости. Рассмотрен расплав свинца при температурах меньше критической (Т < 0.5Тср) и больших отрицательных давлениях и оценена граница его устойчивости.

Публикации автора по теме диссертации

[1] Норман Г. Э., Стегайлов В. В Стохастические свойства молекулярпо-динамической ленард-джонсовской системы в равновесном и неравновесном состояниях // ЖЭТФ. 2001. Т. 119 С. 1011-1020.

(2| Norman G Е., Stegailov V V Stochastic and dynamic properties of molecular dynamics systems Simple liquids, plasma and electrolytes, polymers // Computer Physics Communications 2002 Vol 147. Pp. 678-683

[3] Stegailov V. V. Determinism and Chaos in Decay of Metastable States. Randomness and Predictability of Lifetime // Lecture Notes in Computer Science series (LNCS 2331) / Ed. by P M. A. Sloot, et, al. Springer-Verlag, 2002. Pp. 1147-1153.

[4] Morozov I. V., Norman G E., Stegailov V V. Dynamic and Stochastic Properties of Molecular Dynamics Systems From Simple Liquids to Enzymes // ibid Pp. 1137-1146.

[5] Нормаи Г Э, Стегайлов В. В Гомогенная нуклеация в перегретом кристалле. Молекулярно-дипамический расчет // ДАН. 2002 Т. 386. С. 328-332.

[6] Krivoguz М N, Norman G Е., Stegailov V V., Valuev A A. Superheating of solid metal prior to electric explosion of wires at fast energy deposition ¡/J. Phys. A. 2003. Vol. 36. Pp 6041-6048

[7] Norman G E, Stegailov V. V., Valuev A A. Nanosecond electric explosion of wires: from solid superheating to nonideal plasma formation // Contrib. Plasma Phys. 2003. Vol. 43. Pp. 384-389

[8] Norman G E, Stegailov V. V. Simulation of ideal crystal superheating and decay // Molecular simulation. 2004. Vol. 30. Pp. 397-406.

[9] Kuksm A Yu , Morozov I V, Norman G. E., Stegailov V. V. Standard of molecular dynamics modeling and simulation of relaxation m dense media // Lecture Notes in Computer Science series (LNCS 3039) / Ed by M Bubak, et al. Springer-Verlag, 2004. Pp. 596-603.

[10J Stegailov V. V Homogeneous and heterogeneous mechanisms of superheated solid melting and decay // Computer Physics Communications. 2005. Vol. 169 Pp. 247-250.

[11J Бажиров T T, Норман Г Э, Стерайлов В В Кавитация и область устойчивости жидкого свипца при отрицательных давлениях исследование методом молекулярной динамики // ДАН. 2005. Т. 405, JV« 3.

[12] Stegailov V. V. Optimization of neighbour list techniques and analysis of effects of round-off errors in molecular-dynamics calculations // Beiträge zum Wissenschaftlichen Rechnen Ergebnisse des Gaststudentenprogramms 2002 des John von Neuinann-Tnstituts für Computing / Ed. by R Esser Forschungszentrum Jülich: Zentralinstitut für Angewandte Mathematik, 2002 Pp. 73-86

[13] Стегайлов В В Моделирование процесса нуклеации в перегретом металле при высоких давлениях // Физика экстремальных состояний вещества 2002. Черноголовка: ИПХФ РАН, 2002. С. 46-48.

[14] Стегайлов В В. Молекулярно-дипамическое моделирование плавления в условиях импульсного нагрева // Физика экстремальных состояний вещества 2003 Черноголовка' ИПХФ РАН, 2003. С. 132-133.

[15] Стегайлов В. В. Скорость нуклеации в перегретой кристаллической меди // Физика экстремальных состояний вещества 2005 Черноголовка ИПХФ РАН, 2005 С 164-166.

[16J Бажиров Т. Т., Стегайлов В. В Кавитация в жидком свинце при отрицательных давлениях. Исследование методом молекулярной динамики // там же С 166-168.

Заказ № 2125 Подписано в печать 09.11.2005 Тираж 90 экз. Усл. п л 0,84

, . > ООО "Цифровичок", тел. (095) 797-75-76; (095) 778-22-20 у) www.cfr.ru ; е-тай: mfo@cfr.ru

-7

X

3

200GA

-(£2,5*

Ц- 162Í

Введение

1 Обзор литературы

1.1 Метод молекулярной динамики.

1.2 Равновесие фаз.

1.3 Граница устойчивости фазы

1.4 Кинетика фазовых переходов первого рода.

1.5 Крупномасштабное атомистическое моделирование фазовых превращений при внешних воздействиях

2 Гомогенный распад кристалла

2.1 Термодинамика метастабильных состояний.

2.2 Статистическое описание.

2.3 Метод расчета частоты гомогенной нуклеации.

2.4 Результаты.

2.5 Результаты главы

3 Стохастические свойства метода МД

3.1 Классификация приближений при решении МД задачи

3.2 Неустойчивость траекторий: времена вычислительной и динамической памяти.

3.3 К-энтропия и время динамической памяти.

3.4 Зависимость Ktdm от At и (АЕ2).

3.5 Выбор точности численного интегрирования.

3.6 Физический смысл и роль времени динамической памяти

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

4 Распад и плавление перегретой кристаллической меди

4.1 Гомогенная нуклеация расплава в объеме.

4.2 Температурная зависимость модулей упругости и условия устойчивости.

4.3 Плавление с открытой поверхности

4.4 Результаты главы

5 Кавитации в жидком РЬ при отрицательных давлениях

5.1 Модель и метод расчета.

5.2 Граница устойчивости метастабильной жидкой фазы

5.3 Частота кавитации.

5.4 Обсуждение результатов.

5.5 Результаты главы

Диссертация посвящена разработке метода теоретического исследования динамики и кинетики фазовых превращений в метастабильных кристаллах и жидкостях на молекулярном уровне. Развитый подход применен к неравновесным релаксационным процессам плавления перегретого кристалла и кавитации в растянутой жидкости. Исследованы стохастические свойства метода молекулярной динамики (МД), существенные для рассмотренных задач.

Актуальность работы. Развитие экспериментальной техники, освоение нано-, пико- и фемтосекундных диапазонов открывает возможности получения метастабильных состояний конденсированных веществ со значительными степенями перегрева и/или растяжения по сравнению с равновесным состоянием. Упомянем сильные ударные волны, нагрев фемтосе-кундным лазером, наносекундный электровзрыв проводников и др. Такое развитие науки ставит перед теорией новые задачи, например, о перегреве твердого тела с открытой поверхностью, который ранее в статистической физике считался, вообще говоря, невозможным. Требует более внимательного рассмотрения и классическая теория нуклеации (КТН), надежность которой может оказаться недостаточной, например, при импульсном растяжении жидкостей. Возникающие новые прикладные задачи также стимулируют развитие теории устойчивости метастабильных состояний.

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

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

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

Цель работы состоит в 1) разработке метода МД исследования процессов распада метастабильных состояний с учетом результатов анализа стохастических свойств метода МД, 2) апробации разработанных подходов для расчета кинетики нуклеации в метастабильных кристаллах и жидкостях как для модельных систем, так и для металлов, 3) сравнении полученных результатов с имеющимися теоретическими подходами для описания кинетики фазовых превращений и 4) изучении фазовой диаграммы в области метастабильных конденсированных состояний.

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

Проанализированы стохастические свойства метода молекулярной динамики. Показано, что метод МД является методом, который: i) сохраняет ньютоновскую динамику на времена меньше, чем время динамической памяти, и п) производит статистическое усреднение по начальным условиям вдоль МД траектории. Установлено, что время динамической памяти растет лишь логарифмически при уменьшении уровня шумов. Получены универсальные зависимости, связывающие энтропию Крылова-Колмогорова (максимальный показатель Ляпунова), время динамической памяти и уровень шумов в системе. Введено понятие времени вычислительной памяти — предельного горизонта предсказуемости решения конечно-разностной схемы.

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

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

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

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

Положения, выносимые на защиту.

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] и тезисы российских и международных конференций.

1. Обзор литературы

6.1. Основные результаты и выводы работы

Развиты теоретические основы исследования динамики и кинетики фазовых превращений в метастабильных кристаллах и жидкостях методом

МД.

• Разработан метод расчета частоты гомогенной нуклеации/кавитации, основанный на формировании ансамбля независимых начальных конфигураций и последующем усреднении времени жизни метастабиль-ной фазы по ансамблю МД траекторий (§ 2.3). В качестве апробации метода получена зависимость частоты спонтанной гомогенной нуклеации от степени перегрева кристалла для случаев потенциалов межчастичного взаимодействия мягких сфер и Леннарда-Джонса (рис. 2.5 и 2.6). Проведено МД моделирование распада кристаллической решетки в условиях нагрева с постоянной скоростью. Проверена адекватность описания распада в рамках представления нуклеации как Пуассоновского случайного процесса (рис. 2.4 и 2.9).

• Проанализированы стохастические свойства метода МД и определены его возможности для исследования релаксационных процессов в плотных средах. Получены универсальные логарифмические зависимости, связывающие энтропию Крылова-Колмогорова, время динамической памяти и уровень шумов в системе (формулы (3.14) и (3.15)). Метод МД: i) сохраняет ньютоновскую динамику на времена меньше, чем и ii) производит статистическое усреднение по начальным условиям вдоль МД траектории (§ 3.6). Показано, что в силу принципиальной ограниченности величины ^ возможно лишь статистическое описание распада метастабильного состояния (§ 3.7).

• Рассмотрены гомогенная нуклеация и поверхностное плавление в перегретых кристаллах и кавитация в расплаве при отрицательных давлениях. Показано, что результаты расчетов кинетики гомогенной нуклеации в перегретой кристаллической меди и кавитации в растянутом расплаве свинца находятся в соответствии с оценками по КТН (§ 4.1 и § 5.4). Получена зависимость скорости фронта плавления от температуры (рис. 4.9).

• Найдена граница устойчивости перегретой кристаллической меди при О < Р < 15 ГПа (рис. 4.5): она определяется возникновением сдвиговой неустойчивости кристалла при повышении температуры. Показано соответствие механического и кинетического критериев устойчивости (§ 4.2). Рассмотрен расплав свинца при температурах меньше критической (Т < 0.5Тер) и больших отрицательных давлениях и оценена граница его устойчивости (рис. 5.1).

6.2. Достоверность результатов

Проверялась достоверность полученных результатов. Одни и те же или сходные по смыслу величины рассчитывались разными способами: а) предельно достижимый перегрев Th кристалла в зависимости от скорости нагрева рассчитывался непосредственно и по формуле, связывающей величину Th с частотой спонтанной гомогенной нуклеации J(T), определенной в другом расчете (§ 2.4, п.2); б) установлено совпадение границ устойчивости, найденных из механического и кинетического критериев (§ 4.2).

Результаты, полученные МД расчетом, сопоставлялись с формулами КТН: а) результаты для разного числа частиц согласовывались в рамках предположений КТН (рис. 2.5 и 2.6); б) о связи Th и J(T) сказано выше; в) сопоставление полученных времен жизни перегретой кристаллической меди с оценкой по КТН, дополненной данными о температурной зависимости А(Т) из экспериментов по кристаллизации переохлажденного расплава меди (§4.1).

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

Достоверность моделей проверялась сопоставлением с экспериментальными данными для меди и свинца, путем сравнения результатов расчетов и оценок в рамках КТН (§ 4.1 и § 5.4). Кроме того, температура плавления для меди рассчитывалась двумя независимыми способами (§ 4.3, п.2).

6. Заключение

1. Норман Г. Э., Стегайлов В. В. Стохастические свойства молекулярно-динамической ленард-джонсовской системы в равновесном и неравновесном состояниях // ЖЭТФ.— 2001. — Т. 119.— С. 1011-1020.

2. Norman G. Е., Stegailov V. V. Stochastic and dynamic properties of molecular dynamics systems: Simple liquids, plasma and electrolytes, polymers // Computer Physics Communications. — 2002.— Vol. 147.— Pp. 678-683.

3. Stegailov V. V. Determinism and Chaos in Decay of Metastable States. Randomness and Predictability of Lifetime // Lecture Notes in Computer Science series (LNCS 2331) / Ed. by P. M. A. Sloot, et al. Springer-Verlag, 2002.- Pp. 1147-1153.

4. Норман Г. Э., Стегайлов В. В. Гомогенная нуклеация в перегретом кристалле. Молекулярно-динамический расчет // ДАН. — 2002. — Т. 386. С. 328-332.

5. Krivoguz М. N., Norman G. Е., Stegailov V. V., Valuev A. A. Superheating of solid metal prior to electric explosion of wires at fast energy deposition // J. Phys. A. 2003. - Vol. 36. - Pp. 6041-6048.

6. Norman G. E., Stegailov V. V., Valuev A. A. Nanosecond electric explosion of wires: from solid superheating to nonideal plasma formation // Contrib. Plasma Phys. 2003. - Vol. 43. - Pp. 384-389.

7. Norman G. E., Stegailov V. V. Simulation of ideal crystal superheating and decay // Molecular simulation. — 2004. — Vol. 30. — Pp. 397-406.

8. Stegailov V. V. Homogeneous and heterogeneous mechanisms of superheated solid melting and decay // Computer Physics Communications. — 2005. Vol. 169. - Pp. 247-250.

9. Бажиров Т. Т., Норман Г. Э., Стегайлов В. В. Кавитация и область устойчивости жидкого свинца при отрицательных давлениях, исследование методом молекулярной динамики // ДАН. — 2005. — Т. 405, № 3.

10. Стегайлов В. В. Моделирование процесса нуклеации в перегретом металле при высоких давлениях // Физика экстремальных состояний вещества — 2002. Черноголовка: ИПХФ РАН, 2002. - С. 46-48.

11. Стегайлов В. В. Молекулярно-динамическое моделирование плавления в условиях импульсного нагрева // Физика экстремальных состояний вещества 2003. - Черноголовка: ИПХФ РАН, 2003. - С. 132133.

12. Стегайлов В. В. Скорость нуклеации в перегретой кристаллической меди // Физика экстремальных состояний вещества — 2005. — Черноголовка: ИПХФ РАН, 2005,- С. 164-166.

13. Бажиров Т. Т., Стегайлов В. В. Кавитация в жидком свинце при отрицательных давлениях. Исследование методом молекулярной динамики // Физика экстремальных состояний вещества — 2005. — Черноголовка: ИПХФ РАН, 2005. С. 166-168.

14. Alder В. J., Wainwright Т. Е. Studies in molecular dynamics, i. general method // J. Chem. Phys. ~ 1959. Vol. 31,- Pp. 459-466.

15. Rahman A. Correlations in the motion of atoms in liquid argon // Physical Review. 1964. - Vol. 136. - Pp. A405-A411.

16. Allen M., Tildesley D. Computer Simulation of Liquids. — Oxford: Clarendon press, 1989. — 385 pp.

17. Валуев А., Норман Г., Подлипну к В. Метод молекулярной динамики: теория и приложения // Математическое моделирование. Физико-химические свойства вещества / Под ред. А. А. Самарского, Н. Н. Калиткина. — М.: Наука, 1989. — С. 5-40.

18. Frenkel D., Smith В. Understanding molecular simulation: from algorithms to applications. — San Diego: Academic press, 1996. — 443 pp.

19. Daw M. S., Baskes M. I. Embedded-atom metod: Derivation and application to impurities, surfaces and other defects in metals // Phys. Rev. B. 1984. - Vol. 29. - Pp. 6443-6453.

20. Finnis M. W., Sinclair J. E. A simple empirical N-body potential for transition metals // Phil. Mag. A. — 1984. Vol. 50. - Pp. 45-55.

21. Stillinger F. H., Weber T. A. Computer simulation of local order in condensed phases of silicon // Phys. Rev. B. 1985. - Vol. 31. - Pp. 52625271.

22. Ercolessi F., Parrinello M., Tosatti E. Simulation of gold in the glue model // Phil Mag. A. 1988. - Vol. 58. - Pp. 213-226.

23. Tersoff J. New empirical approach to the structure and energy of covalent systems 11 Phis. Rev. B. 1988. - Vol. 37. - Pp. 6991-7000.

24. Lim H. S., Ong С. K., Ercolessi F. Stability of face-centered cubic and icosahedral lead clusters // Surface Science. — 1992,— Vol. 269/270. — Pp. 1109-1115.

25. Mishin Y., Mehl M. JPapaconstantopoulos D. A., Voter A. F., Kress J. D. Structural stability and lattice defects in copper: Ab initio, tight-binding, and embedded-atom calculations // Phys. Rev. B. — 2001,-Vol. 63.- P. 224106.

26. Mendelev M. I., Han S., Srolovitz D. J., Ackland G. J., Sun D. Y., As-ta M. Development of new interatomic potentials appropriate for crystalline and liquid iron // Phil. Mag. 2003. - Vol. 83. - Pp. 3977—3994.

27. Белащенко Д. К. Семейства межмолекулярных потенциалов, приводящих к тождественным структурам некристаллических тел в методе молекулярной динамики // Журнал физической химии. — 2004. — Т. 78.-С. 1621-1628.

28. Ландау JI. Д., Лифшиц Е. М. Теоретическая физика. — М.: Физмат-лит, 2001. — Т. V. Статистическая физика. — 616 с.

29. Левашов П. Р. Уравнения состояния жидкой фазы при высоких давлениях и температурах // Препринт ОИВТ РАН № 1-446. — 2000. — 29 с.

30. Hoover W. G. Canonical dynamics: Equilibrium phase-space distributions 11 Phys. Rev. A. 1985. - Vol. 31. - Pp. 1695-1697.

31. Berendsen H. J. C., Postma J. P. M., van Gunsteren W. F., DiNola A., Haak J. R. Molecular dynamics with coupling to an external bath // «7. Chem. Phys. 1984. - Vol. 81. - Pp. 3684-3690.

32. Robbins M. О., Grest G. S., Kremer K. Effect of finite system size on thermal fluctuations: Implications for melting // Phys. Rev. B. — 1990. — Vol. 42. Pp. 5579-5585.

33. Румер Ю. Б., Рывкин M. Ш. Термодинамика, статистическая физика и кинетика. — Новосибирск: Изд-во Новосиб. ун-та, 2000. — 608 с.

34. Broughton J. Q., Gilmer G. H. Molecular dynamics investigation of the crystal-fluid interface. I. Bulk properties // J. Chem. Phys. — 1983. — Vol. 79. Pp. 5095-5104.

35. Lutsko J. F., Wolf D., Phillpot S. R., Yip S. Molecular dynamics study of lattice-defect-nucleated melting in metals using an embedded-atom-method potential // Phys. Rev. B. 1989. - Vol. 40. - Pp. 2841-2855.

36. Morris J. R., Wang C. Z., Но К. M., Chan С. T. Melting line of aluminum from simulations of coexisting phases // Phys. Rev. B. — 1994. — Vol. 49.- Pp. 3109-3115.

37. Belonoshko А. В., Ahuja R., Johansson B. Quasi-ab initio molecular dynamic study of Fe melting // Phys. Rev. Lett. — 2000.— Vol. 84.— Pp. 3638-3641.

38. Baidakov V. G., Protsenko S. P., Chernykh G. G., Boltachev G. Sh. Statistical substantiation of the van der Waals theory of inhomogeneous fluids 11 Phys. Rev. E. 2002. - Vol. 65. - P. 41601.

39. Скрипов В. П., Коверда В. П. Спонтанная кристаллизация переохлажденных жидкостей. — М.: Наука, 1984. — 232 с.

40. Борн М., Хуан К. Динамическая теория кристаллических решеток. — М.: Изд. ин. лит., 1958.-488 с.

41. Байдаков В. Г., Галашев А. Е., Скрипов В. П. Устойчивость перегретого кристалла в молекулярно-динамической модели аргона // Физ. те. тела. ~ 1980. Т. 22. - С. 2682-2686.

42. Wang J., Li J., Yip S., Phillpot S., Wolf D. Mechanical instabilities of homogeneous crystals // Phys. Rev. В. 1995.- Vol. 52.- Pp. 1262712635.

43. Wang J., Li J., Yip S., Wolf D., Phillpot S. R. Unifying two criteria of born: Elastic instability and melting of homogeneous crystals // Physica A. 1997. - Vol. 240. - Pp. 396-403.

44. Jin Z. H., Gumbsch P., Lu К., Ma E. Melting mechanisms at the limit of superheating 11 Phys. Rev. Lett. 2001. - Vol. 87. - P. 55703.

45. Зельдович Я. Б. К теории образования новой фазы. Кавитация // ЖЭТФ. 1942. - Т. 12, № 11/12. - С. 525-538.

46. Скрипов В. П. Метастабильная жидкость. — М.: Наука, 1972. — 312 с.

47. Motorin V. I., Musher S. L. Kinetics of the volume melting. Nucleation and superheating of metals // J. Chem. Phys.— 1984.— Vol. 81.— Pp. 465-469.

48. Lu K., Li Y. Homogeneous nucleation catastrophe as a kinetic stability limit for superheated crystal // Phys. Rev. Lett. — 1998. — Vol. 80. — Pp. 4474-4477.

49. Iwamatsu M., Horii K. A one order parameter theory of crystal-melt nucleation of spherical clusters // J. Phys. Soc. Japan. — 1996. — Vol. 65. — Pp. 2311-2316.

50. Iwamatsu M. Homogeneous nucleation for superheated crystal // J. Phys.: Condens. Matter. 1999. - Vol. 11. - Pp. L1-L5.

51. Ma D., Li Y. Heterogeneous nucleation catastrophe on dislocations in superheated crystals // J. Phys.: Condens. Matter. — 2000. — Vol. 12.— Pp. 9123-9128.

52. Kinjo Т., Matsumoto M. Cavitation processes and negative pressure // Fluid Phase Equilibria. 1998. - Vol. 144. - Pp. 343-350.

53. Жуховицкий Д. И. Поверхностное натяжение границы раздела пар-жидкость с конечной кривизной // Коллоидный журнал. — 2003. — Т. 65. С. 480-494.

54. Luo S.-N., Ahrens Т. J., Qagin Т., Strachan A., Ill W. A. G., Swift D. C. Maximum superheating and undercooling: systematics, molecular dynamics simulations, and dynamic experiments // Phys. Rev. В.— 2003.— Vol. 68.-P. 134206.

55. Zhakhovskii V. V., Zybin S. V., Nishihara K., Anisimov S. I. Orientation dependence of shock structure with melting in L-J crystal from molecular dynamics // Prog. Theor. Phys. Suppl 2000,- Vol. 138.- Pp. 223228.

56. Анисимов С. И., Жаховский В. В., Иногамов Н. А., Нишихара К., Опарин А. М., Петров Ю. В. Разрушение твердой пленки в результате действия ультракороткого лазерного импульса // Письма в ЖЭТФ. 2003. - Т. 77. - С. 731-736.

57. Tanguy D., Mareschal М., Lomdahl P. S., Germann Т. С., Holian В. L., Ravelo R. Dislocation nucleation induced by a shock wave in a perfect crystal: Molecular dynamics simulations and elastic calculations // Phys. Rev. B. 2003. - Vol. 68. - P. 144111.

58. Abraham F. F. How fast can cracks move? A research adventure in materials failure using millions of atoms and big computers // Advances in Physics. 2003. - Vol. 52. - Pp. 727—790.

59. Ivanov D. S., Zhigilei L. V. Combined atomistic-continuum modeling of short-pulse laser melting and disintegration of metal films // Phys. Rev. B. 2003. - Vol. 68. - P. 64114.

60. Wood G. R., , Walton A. G. Homogeneous nucleation kinetics of ice from water // J. Appl. Phys. 1970. - Vol. 41. - Pp. 3027-3036.

61. Hoover W. G., Gray S. G., Johnson K. W. Thermodynamic properties of the fluid and solid phases for inverse power potentials // J. Chem. Phys. — 1971. Vol. 55. - Pp. 1128-1136.

62. Кривогуз M. H., Норман Г. Э. Спинодаль перегретого твердого металла // ДАН. 2001. - Vol. 379. - Pp. 177-180.

63. Hoover W., Ross M., Johnson K. W., Henderson D., Barker J. A., Brown В. C. Soft-sphere equation of state // J. Chem. Phys. — 1970. — Vol. 52,- Pp. 4931-4941.

64. Hardy W. H., Crawford R. K., Daniels W. B. Experimental determination of the P-T melting curve for argon //J. Chem. Phys. — 1971. — Vol. 54. — Pp. 1005-1010.

65. Barroso M. A., Ferreira A. L. Solid-fluid coexistence of the Lennard-Jones system from absolute free energy calculations //J. Chem. Phys. —2002. Vol. 116. - Pp. 7145-7150.

66. Скрипов В. П., Файзуллин М. 3. Фазовые переходы кристалл-жидкость-пар и термодинамическое подобие.— М.: Физматлит,2003.- 160 с.

67. Rethfeld В., Sokolowski-Tinten К., von der Linde D., Anisimov S. I. Ultrafast thermal melting of laser-excited solids by homogeneous nucleation 11 Phys. Rev. B. 2002. - Vol. 65. - P. 92103.

68. Заславский Г. M. Стохастичность динамических систем. — М.: Наука, ГРФМЛ, 1984.- 272 с.

69. Ландау Л. Д., Лифшиц Е. М. Теоретическая физика. — М.: Физмат-лит, 2001. — Т. III. Квантовая механика (нерелятивистская теория). — 808 с.

70. Кравцов Ю. А. Фактические границы гипотезы замкнутости и классические парадоксы кинетической теории // ЖЭТФ. — 1989. — Т. 96. — С. 1661-1665.

71. Кравцов Ю. А. Случайность, детермированность, предсказуемость // УФН. 1989. - Т. 158. - С. 93-122.

72. Кравцов Ю. А. Фундаментальные и практические пределы предсказуемости // Пределы предсказуемости / Под ред. Ю. А. Кравцова.— Москва: ЦентрКом, 1997. С. 170-200.

73. Герценштейн М. Е., Кравцов Ю. А. Ограничения применимости ньютоновского описания движения частиц в газе вследствие спонтанного излучения низкочастотных фотонов // ЖЭТФ, — 2000,— Т. 118.— С. 761-763.

74. Малинецкий Г. Г. Синергетика, предсказуемость и динамический хаос // Пределы предсказуемости / Под ред. Ю. А. Кравцова.— М.: ЦентрКом, 1997. С. 78-139.

75. Климонтович Ю. Л. Статистическая теория открытых систем. — М.: Янус, 1995. 624 с.

76. Hoover W. G. Time Reversibility, Computer Simulation and Chaos. — Singapore: World Scientific, 1999. 280 pp.

77. Hoover W. G. Computational Statistical Mechanics. — Amsterdam: Elsevier, 1991.-314 pp.

78. Posch H. A., Hoover W. G. Lyapunov instability of dense Lennard-Jones fluids 11 Phys. Rev. A. 1998. - Vol. 38. - Pp. 473-482.

79. Hoover W. G., Posch H. A. Second-law irreversibility and phase-space dimensionality loss from time-reversible nonequilibrium steady-state Lya-punov spectra // Phys. Rev. E. 1994. - Vol. 49. - Pp. 1913-1920.

80. Kwon K.-H., Park B.-Y. Lyapunov exponent and the solid-fluid phase transition // J. Chem. Phys. 1997. - Vol. 107.-Pp. 5171-5178.

81. Каклюгин А. С., Норман Г. Э. Hierarchical approach generalization of vitalism and reductionism // Российский химический журнал. — 2000.-Т. 44(3).-С. 7-20.

82. Валуев А. А., Норман Г. Э., Подлипчук В. Ю. Уравнения метода молекулярной динамики // Термодинамика необратимых процессов / Под ред. А. И. Лопушанская. — М.: Наука, 1987. — С. 11-17.

83. Norman G. Е., Podlipchuk V. Yu., Valuev A. A. On the theory of the molecular-dynamics method 11 J. Moscow Phys. Soc. — 1992. — Vol. 2. — Pp. 7-21.

84. Валуев А. А., Норман Г. Э., Подлипчук В. Ю. Энтропия Крылова-Колмогорова неупорядоченных Леннард-Джонсовских систем // Математическое моделирование. — 1990. — Т. 2(5). — С. 3-7.

85. Loskutov A. Chaotic dynamics of chemical systems // Mathematical Methods in Contemporary Chemistry / Ed. by S. I. Kuchanov. — Amsterdam: Gordon and Breach, 1996.— Pp. 181-265.

86. Метод молекулярной динамики в физической химии / Под ред. Ю. К. Товбина. М.: Наука, 1996. - 334 с.

87. Белащенко Д. К. Механизмы диффузии в неупорядоченных системах (компьютерное моделирование) // УФЕ. — 1999.— Т. 169.— С. 361— 384.

88. Morozov I. V., Norman G. E., Valuev A. A. K-entropies of electrons and ions in nonideal plasmas // Contrib. Plasma Phys. — 1999. — Vol. 39. — Pp. 307-311.

89. Morozov I. V., Norman G. E., Valuev A. A. K-entropy (average lyapunov exponent), dynamics and chaos for particle // J. de Physique (France). — 2000. Vol. 10(Pr5). - Pp. 251-254.

90. Morozov I. V., Norman G. E., Valuev A. A. Stochastic properties of strongly coupled plasmas // Phys. Rev. E. 2001. - Vol. 63. - P. 36405.

91. Ueshima Y., Nishihara K., Barnett D. M., Tajima Т., Furukawa H. Particle simulation of Lyapunov exponents in one-component strongly coupled plasmas // Phys. Rev. E. 1997. - Vol. 55,- Pp. 3439-3449.

92. Ueshima Y., Nishihara K., Barnett D. M., Tajima Т., Furukawa H. Relation between Lyapunov exponent and dielectric response function in dilute one component plasmas // Phys. Rev. Lett. — 1997.— Vol. 79.— Pp. 2249-2252.

93. Barnett D. M., Tajima Т., Nishihara K., Ueshima Y., Furukawa H. Lyapunov exponent of a many body system and its transport coefficients // Phys. Rev. Lett. 1996. - Vol. 76. - Pp. 1812-1815.

94. Norman G. Е., Podlipchuk V. Yu., Valuev A. A. Equations of motion and energy conservation in molecular dynamics // Molecular Simulation. — 1993. Vol. 9. - Pp. 417-424.

95. Rowlands G. A numerical algorithm for Hamiltonian systems // J. Computational Physics. 1991. - Vol. 97. - Pp. 235-239.

96. Lopez-Marcos M. A., Sanz-Serna J. M., Diaz J. C. Are Gauss-Legendre methods useful in molecular dynamics? // J. Comput. Appl. Math. —1996. Vol. 67. - Pp. 173-179.

97. Lopez-Marcos M. A., Sanz-Serna J. M., Skeel R. D. Explicit symplectic integrators using Hessian-vector products / / SI AM J. Sci. Comput. —1997. Vol. 18. - Pp. 223-238.

98. Kaklyugin A. S., Norman G. E. Quantum corrections to the classical equations of motion //J. Moscow Phys. Soc. — 1995. — Vol. 5. — Pp. 167— 180.

99. Рудяк В. Я., Харламов Г. ВБелкин А. А. Диффузия наночастиц и макромолекул в плотных газах и жидкостях // ТВ Т.— 2001.— Т. 39.- С. 283-291.

100. Коверда В. П., Скрипов В. П. Спонтанная кристаллизация переохлажденных жидких металлов // ФММ. — 1973. — по. 35. — Pp. 988-992.

101. Ray J. R. Elastic constants and statistical ensembles in molecular dynamics 11 Computer Physics Reports. — 1988. Vol. 8. - Pp. 109—152.

102. Wolf R. J., Mans our K. A., Lee M. W., Ray J. R. Temperature dependence of elastic constants of embedded-atom models of palladium // Phys. Rev. B. 1992. - Vol. 46. - Pp. 8027—8035.

103. Nguyen Т., Ho P. S., Kwok Т., Nitta C., Yip S. Thermal structural disorder and melting at a crystalline interface // Phys. Rev. В. — 1992. — Vol. 46. Pp. 6050-6060.

104. Bilalbegovic G., Lutz H. 0. The onset of a liquid-vapour transition in metallic nanoparticles // Chemical Physics Letters. — 1997. — Vol. 280. — Pp. 59-65.

105. Ternovoi V. Y., Fortov V. E., Kvitov S. V., Nikolaev D. N. Experimental study of lead critical point parameters // Shock Compression of Condensed Matter (Part 1) / Ed. by S. C. Schmidt, W. C. Tao. New York: AIP Press, 1996. - Pp. 5-40.

106. Медин С. А., Орлов Ю. H., Паршиков A. H., Суслин В. M. Моделирование отклика первой стенки камеры и бланкета реактора ИТС на микровзрыв (препринт №41).— М.: ИПМ им. М.В.Келдыша РАН, 2004. 32 с.

107. Иосилевский И. JI., Чигвинцев А. Ю. Спинодальный распад зоны метастабильного плавления в пределе нулевой температуры // Электронный журнал иИсследовано в России" (URL: <http://zhurnal.ape.relarn.ru/articles/2003/003.pdf>). 2003. - С. 2134.

108. Бабичев А., Бабушкина Н. А., Братковский А. М., и др. Физические величины: Справочник / Под ред. И. С. Григорьева, Е. 3. Мейлихо-ва. — М.: Энергатомиздат, 1991.— 1232 с.