Авторефераты по всем темам  >>  Авторефераты по физике

На правах рукописи

УДК 538.9 Стегайлов

Владимир Владимирович Классические и квантовые атомистические модели отклика конденсированных сред на интенсивные энергетические воздействия

Специальность 01.04.07 физика конденсированного состояния

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

Москва 2012

Работа выполнена в Федеральном государственном бюджетном учреждении науки УОбъединенный институт высоких температур РАНФ на кафедре физики высокотемпературных процессов Московского физико-технического института (Государственного университета).

Научный консультант: доктор физико-математических наук, заведующий отделом, профессор Генри Эдгарович Норман (Объединенный институт высоких температур РАН)

Официальные оппоненты: доктор физико-математических наук, ведущий научный сотрудник, Наиль Алимович Иногамов (Институт теоретической физики им. Л.Д.Ландау РАН) доктор физико-математических наук, заместитель директора по науке, Валентин Николаевич Рыжов (Институт физики высоких давлений им. Л.Ф.Верещагина РАН) доктор физико-математических наук, главный научный сотрудник, профессор Вячеслав Михайлович Чернов (ОАО Высокотехнологический научно-исследовательский институт неорганических материалов им. академика А.А.Бочвара )

Ведущая организация: Институт прикладной математики им. М.В.Келдыша РАН

Защита состоится 2012 г. в часов на заседании диссертационного совета Д.212.156.06 при Московском физико-техническом институте (Государственном университете) по адресу: 117393, г. Москва, ул. Профсоюзная, д. 84/32, корпус В-2.

Отзывы направлять по адресу: 141700, г. Долгопрудный Московской обл., Институтский переулок, д. 9., МФТИ, Диссертационный совет Д 212.156.06.

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

Автореферат разослан 2012 г.

Ученый секретарь диссертационного совета кандидат технических наук Н. П. Чубинский 1

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

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

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

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

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

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

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

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

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

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

Цель работы состоит: 1) анализе теоретических основ метода молекулярной динамики, 2) разработке межатомных потенциалов взаимодействия, 3) моделировании различных случаев нуклеации новой фазы при фазовых переходах первого рода, 4) построении моделей пластической деформации на основе динамики одиночных дислокаций, 5) построении моделей ударноволнового разрушения твердых тел, 6) радиационных повреждений и 7) абляции металлов.

Научная новизна работы. Проанализированы стохастические ограничения на расчет динамических траекторий в методе МД. Показано, что метод МД является методом, который: i) сохраняет ньютоновскую динамику на времена меньше, чем время динамической памяти, и ii) производит статистическое усреднение по начальным условиям вдоль МД траектории.

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

Развита методика создания эффективного многочастичного потенциала межатомного взаимодействия исключительно на основе ab initio расчетов.

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

Проведено исследование кинетики объемного и поверхностного плавления, установлен столкновительно-лимитированный тип кинетики движения фронта плавления, выявлен эффект предплавления металлов в контакте с неупорядоченной средой.

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

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

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

Ab initio расчеты в рамках теории функционала электронной плотности использованы для анализа устойчивости твердых металлов (Au) и диэлектриков (LiF) в состоянии разогретого плотного вещества после импульсного энерговклада в электронную подсистему. Развита двухтемпературная модель абляции металлов под действием субпикосекундных импульсов. Обнаружен механизм абляции, связанный с релаксацией электронного давления.

Практическая ценность работы.

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

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

1. Стохастические ограничения предсказуемости динамических корреляций в методе молекулярной динамики. Время динамической и вычислительной памяти. Усредненно-статистический смысл молекулярно-динамических траекторий.

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

3. Температурная зависимость подвижности краевой дислокации в монокристаллическом алюминии.

4. Столкновительно-лимитированный характер кинетики движения фронта плавления в металлах. Эффект предплавления металлов в контакте с внешней средой-медиатором.

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

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

7. Метод создания межатомных потенциалов для металлов на основе первопринципных расчетов. Межатомный потенциал для системы Mo-Xe.

Межатомный потенциал для золота, зависящий от электронной температуры.

8. Сильное влияние электронной температуры на устойчивость кристаллов LiF и Au в состоянии разогретого плотного вещества.

9. Различие двух механизмов абляции золота при воздействии субпикосекундных импульсов: механизма, связанный с релаксацией электронного давления в случае двухтемпературного металла с горячей электронной подсистемой, и механизма, связанного с распространением волны разгрузки в равновесном металле.

Апробация работы. Основные результаты диссертационной работы были доложены на конференциях УУравнения состояния веществаФ (п. Эльбрус, 2002, 2004 и 2010), УВоздействие интенсивных потоков энергии на веществоФ (п. Эльбрус, 2003, 2005 и 2011), УNucleation theory and applicationsФ (Дубна, 2003-2008), УFoundations of molecular modeling and simulationФ (США, Keystone, 2003), УФазовые превращения при высоких давленияхФ (Черноголовка, 2004), Computational Physics (Германия, Aachen, 2001; США, San Diego, 2002;

Италия, Genova, 2004; Gyeongju, Республика Корея, 2006), 15Цth Symposium on thermophysical properties (США, Bolder, 2003), научно-координационных сессиях УИсследования неидеальной плазмыФ (Президиум РАН, Москва, 2001, 2004, 2010), симпозиумах УПроблемы физики ультракоротких процессов в сильнонеравновесных средахФ (Новый Афон, 2003-2011), Конференциях Американского физического общества УShock compression of condensed matterФ (Гавайи, США, 2007; Нэшвил, США, 2009; Чикаго, США, 2011), Международных семинарах УComplex Systems of Charged Particles and their Interaction with Electromagnetic RadiationФ (Москва, 2008-2010), Международном семинаре УNew Models and Hydrocodes for Shock Wave Processes in Condensed MatterФ (Эшторил, Испания, 2008), Международных форумах по нанотехнологиям (Москва, 2008, 2009, 2010), Международной конференции УStatistical Physics: Modern Trends and Applications" (StatPhys2009) (Львов, Украина, 2009), Российско-японском семинаре УMolecular Simulation Studies in Materials and Biological SciencesФ (Дубна, 2010), Конференциях УФизика сильно сжатого веществаФ (Троицк, 2009, 2010).

Публикации. По материалам диссертации опубликовано 27 работ в реферируемых научных изданиях.

Структура и объем диссертации. Диссертация состоит из введения и семи глав, изложена на 217 страницах, включая 96 рисунков и 381 наименований цитируемой литературы.

4 3 1 2 0 3 4 0 2 4 0. 0.3 0.02 10- 10-1 t'm 10-t Г0 Г01(a) (б) Рис. 1: Слева: Иллюстрация начального участка траектории, рассчитываемой в ММД из начального состояния системы 0. Сплошной кривой показано точное решение соответствующей задачи Коши. Пунктирными кривыми показаны Ньютоновские траектории МД системы, на которые попадает решение конечно-разностной аппроксимации.

Справа: Нормированные усредненные разбегания координат r2(t) и скоростей v2(t) на двух траекториях, рассчитанных из тождественных начальных условий с шагами t=0.001 и t =0.0001 (N = 64, = 0.5, T = 0.44, потенциал Леннарда-Джонса).

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

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

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

Проанализировано соотношение между динамическими и стохастическими свойствами траекторий частиц, обусловленное Ляпуновской неустойчивостью уравнений Ньютона и ошибками схемы и округления численного интегрирования. Рассмотрено понятие времени динамической памяти (или предсказуемости) td для траекторий частиц. Метод МД приближенно сохраняет m 1t'm d Ktm t = 0.0E / E t' / t 1x10-25 1x10-20 1x10-15 1x10-10 1x10-0 0.1 0.2 0.3 0.4 0.Рис. 2: Слева: Зависимость tm от t /t для различных значений t: 1 - t = 0.01, 2 - 0.005, 3 - 0.001 (крайняя левая точка на линии 1 соответствует t = 0.00001); (N = 64, = 0.5, T = 0.44, потенциал Леннарда-Джонса).

Справа: Предполагаемое ограничение предельных значений Ktd точностью представлеm ния вещественных чисел в компьютере: 1, 2 и 3 соответствуют 7, 16 и 31 десятичным знакам в машинном представлении действительных чисел. Зависимость безразмерного произведения Ktd от относительной флуктуации полной энергии системы: кружки и m кресты - Леннард-Джонсовская система и неидеальная плазма соответственно (неявная схема Эйлера), ромбы - Леннард-Джонсовская система (схема Рунге-Кутты 4-го порядка точности). Пунктир 4 соответствует (1).

Ньютоновскую динамику только на временах, меньших времени td, а даm лее траектории расходятся. Образуется пучок траекторий на Ньютоновской УножкеФ вместо единственного решения задачи Коши (рис. 1 и рис. 2).

Траектории частиц, вычисляемые в ММД, соответствуют уравнениям движения, которые отличаются от Ньютоновских малыми слагаемыми. Эти слагаемые имеют стохастический характер и, таким образом, системы частиц, рассматриваемые в ММД, не являются Гамильтоновыми. Вместо микроканонического в ММД возникает ансамбль NV (E E), где E флуктуация полной энергии E вдоль МД траектории. Такой ансамбль ранее вводился в эргодической теории.

Результаты метода МД содержат все базовые понятия статистической физики: необратимость, переход к статистическому описанию, а для равновесных систем Гиббсовское распределение и равенство средних по времени средним по пространству. Из анализа того, как эти результаты получены следуют важные выводы. Невозможно решить систему уравнений Ньютона в течение времен, требуемых для усреднения в статистической физике, Рис. 3: Парные потенциалы: (a) сплошная красная линия Mo-Mo; черный пунктир Mo-Xe; синий короткий пунктир Xe-Xe. Эффективная электронная плотность (b) и функции погружения (c): сплошная красная линия Mo; синий пунктир Xe.

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

Ktd = - ln( E2 ) + const. (1) m Увеличение времени динамической памяти путем повышения точности схемы ограничено конечной точность представления вещественных чисел в расчетах (т.е. ошибками округления, которые определяют время вычислительной памяти, см. рис. 2).

Из найденных решений фундаментальных проблем следуют стандарты требований, которым должны удовлетворять МД расчеты для конкретных задач как при исследовании систем, выведенных на равновесие, так и при изучении релаксации. Выбор числа частиц N определяется физическими факторами размерами системы, характерными для пространственных и временных масштабов задачи, которую предполагается исследовать. Выбор N, в свою очередь, определяет оптимальное число вычислительных ядер (процессоров), превышение которого нецелесообразно при параллельных вычислениях. Точность усреднения в равновесной МД не хуже, чем (td /t0)1/2, m где t0 длительность МД траектории. Для релаксационных процессов усреднение проводится по ансамблю статистически независимых микроскопических начальных состояний, эквивалентных друг другу макроскопически.

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

Идея метода согласования по силе (force-matching method) заключается в построении потенциала из ab initio данных и, в первую очередь, по вычисленным силам, действующим между атомами. Развития межатомных потенциалов идет в направлении полного перехода на данные первопринципных расРис. 4: Слева: Кривая теплового расширения: звездочки - экспериментальные данные из (Edwards et al., 1951); кривые - зависимости полученные в молекулярно-динамическом моделировании с использованием различных потенциалов.

Справа: Изотерма молибдена при комнатной температуре: экспериментальные данные и результаты расчетов с различными потенциалами.

четов. В данной работе для квантовых расчетов использовался программный код VASP (Kresse and Furthmuller, 1996), основанный на методе функционала плотности (DFT) и базисе плоских волн. Квантовые расчеты, использовались как входные данные в методе force-matching. Для подгонки потенциальных параметров использовался программный код PotFit (Brommer and Gahler, 2007), где для поиска потенциала применяется метод блуждания броуновской частицы с понижающейся температурой в градиентном поле многомерного пространства параметров.

Первая стадия разработки потенциала заключалась в создании систем с относительно небольшим числом атомов (около 50 - 300), находящихся в различных фазах и состояниях (для системы молибден-ксенон использовалась 81 конфигурация). После создания набора конфигураций следовала стадия квантовых расчетов. Для каждой конфигурации проводился расчет с помощью метода функционала электронной плотности. Полученные значения являлись эталонными значениями в алгоритме force-matching.

Третья стадия заключалась в поиске многочастичного потенциала погруженного атома (EAM-потенциала), который давал бы значения энергий, напряжений и сил максимально близких к эталонным значениям. В случае EAM-потенциала суммарная энергия системы определяется как U = (rij) + F(j). (2) i

C <111> D <111> D <110> вакансия DFT (Han at el., 2002) 7.34 7.34 7.51 DFT (Dudarev at el., 2006) 7.419 7.417 7.581 DFT (Finkenstadt, 2006) - 7.70 7.85 МД (Dudarev et al., 2007) 7.432 7.406 7.582 МД (Finnis-Sinclair, 1987) 7.148 7.148 6.944 МД (Zhang Wang Xu, 2009) 5.54 5.48 5.83 МД (новый потенциал) 6.42 6.43 6.67 2.функции для Mo-Mo, Mo-Xe и Xe-Xe взаимодействий). Второе слагаемое учитывает многочастичное взаимодействие, поскольку функция УпогруженияФ F зависит от эффективной Уэлектронной плотностиФ j в месте расположения j-го атома, наводимой всеми остальными атомами системы:

j = (rij), (3) i =j где (rij) представляет собой функцию Уэлектронной плотностиФ наводимой атомом i (типа ). Для каждого из типа атомов существуют свои функции F и (rij). Все 7 независимых функций задаются кубическими сплайнами.

Ошибка подгонки представляется как Z-функция состоящая из двух частей:

Z = ZF + ZC, (4) где Na EAM DF (Fi,x - Fi,x T )j j ZF = Wi DF T (5) (Fi,x )2 + j i=1 j=и Nc (AEAM - ADF T )i ZC = Wi i. (6) (ADF T )2 + i i=В этих выражениях введены следующие обозначения: Na общее число всех атомов во всех конфигурациях, Nc общее число энергий систем и компонент тензора напряжений, Fi,x xj компонента силы действующей на i j атом, Wi статистический вес i конфигурации, Ai энергия или компонента тензора напряжений в i конфигурации, малая величина предотвращающая переоценку конфигурации с маленькими эталонными значениями, EAM и DFT значения рассчитанные для данной конфигурации с помощью EAM-потенциала и эталонные значения рассчитанные квантовыми методами соответственно.

Рис. 5: Слева: Основные типы дефекта межузельного атома в оцк-кристалле (например, молибдене): а D <111>; б D <110>; в D <100>; г C <111> (crowdion); д T (tetrahedral); е O (octahedral). По данным ab-initio расчетов, дефекты типа а и г должны обладать наименьшей энергией образования (быть наиболее стабильными).

Справа: Энергетическая иерархия SIA-дефектов. Энергия образования отсчитана от энергии образования дефекта D<111>. Приведены результаты двух различных квантовых расчетов (Han et al., 2002; Nguyen-Manh et al., 2006) и результаты моделирования с двумя различными EAM-потенциалами.

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

В третей главе рассмотрена кинетика элементарных процессов при фазовых переходах первого рода на примере плавления и кавитации. Представлена многомасштабная модель Унуклеация и ростФ для описания откола в жидкости.

Объемное плавление. В большинстве случаев плавление твердого тела начинается на неоднородностях кристаллической структуры: дефектах, межзеренных границах, открытой поверхности. Если при некоторых особых условиях роль гетерогенных механизмов плавления практически устраняется (например, в объеме монокристаллов), плавление происходит в результате гомогенной нуклеации, т.е. случайного флуктуационного образования зародыM P,GPa mpmpS T,K 800 1200 1600 2000 24M P,GPa mmmpmp1 ppS T,K 800 1200 1600 2000 24Рис. 6: Иллюстрация процесса изохорического нагрева в координатах давление - температура. M - линия плавления меди, хорошо воспроизводящаяся EAM моделью.

mmS - граница устойчивости кристалла. Показан срез атомной структуры в расчетной ячейке (более светлые атомы имеют менее упорядоченное локальное окружение): p1 и p2 - поликристаллическое твердое тело до начала до начала зернограничного плавления и после 250 пс, m1 - сильно перегретая решетка, m2 - рост зародыша, образовавшегося в результате гомогенной нуклеации. Размер системы 0.5 миллиона атомов.

p1 pшей расплава. Гетерогенное плавление происходит практически безбарьерно (и является основным ограничением на перегрев кристаллической фазы в эксперименте), для начала гомогенного зародышеобразования требуются достаточно большие степени перегрева, соответствующие высоким энергиям активации. Такие перегревы могут реализовываться в условиях импульсных энерговкладов большой мощности за времена недостаточные для существенного развития гетерогенного плавления.

Пример МД системы для моделирования плавления в объеме представляет собой кубическую расчетную ячейку в трехмерных периодических граничных условиях. Используется EAM потенциал для меди (Mishin et al., 2001).

Исходная структура соответствует монокристаллу или поликристаллической структуре, выведенной на равновесие при нулевом давлении. Поликристаллическая структура создается разделение пространства расчетной ячейки на многогранники Вороного и последующим их заполнением г.ц.к. структурами со случайной ориентацией (возникающие сильно накладывающиеся атомы на границах зерен удаляются до начала МД расчета). Система нагревается с использованием Ланжевеновского термостата со скоростью 2.3 К/пс (см.

рисунок 6). Два указанных случая рассматриваются в качестве предельных случаев, описывающих поликристаллические твердые тела с наноразмерными зернами и с крупными (например, микронными) зернами. Конечная скорость распространения фронта плавления с межзеренных границ может привести к перегреву объема зерен, степень которого ограничивается гомогенной vfront /vth 0.Рис. 7: Зависимость vfront(T ) в безразмерном виде vfront/vth и T/T. Одновременно приведены результаты для (100), (110) и (111) ориентаций. Результаты для меди (0.LJ - данная работа, 2 - (Lutsko et al., 1989)) Cu получены подобным образом. Пунктирная Al Fe линия показывает единичный наклон.

Cu 0 0.05 0.1 0.15 0.T/T нуклеацией вблизи границы устойчивости. Этот эффект исчезает когда размер зерен уменьшается. Максимально достижимые степени перегрева могут также зависеть от электронной теплопроводности, которая определяет подвод тепла в область фронта плавления (что не учитывается в данной модели).

Поверхностное плавление. В данном главе приводятся результаты по зависимости скорости фронта плавления от температуры для Леннард-Джонсовского кристалла, кристаллического алюминия и железа, которые моделируются на основе EAM потенциалов. Для описания столкновительно-лимитированной кинетики роста при кристаллизации была предложена теория Броутона-Гилмера-Джексона. В переделе малых степеней переохлаждения T скорость кристаллизации зависит от температуры как vfront vthT/T, где vth = (3kBT/m)1/2 - тепловая скорость и T = T - Tm (Jackson, 2004).

Результаты расчетов показывают, что БГД-теория кристаллизации может быть перенесена на случай плавления. На рисунке 7 результаты, полученные для различных систем представлены в безразмерных переменных vfront/vth и T/T. Зависимости имеют почти одинаковый наклон, близкий к единице (с небольшими отличиями для Al). Это свидетельствует, что кинетика плавления в моноатомных твердых телах является столкновительнолимиторованным процессом и его скорость может быть оценена как vthT/T вплоть до значений перегрева T/T = 0.2.

Поверхностное предплавление металлов при высоких давлениях. Были выполнены МД расчеты по изучению процесса поверхностного плавления железа, молибдена и алюминия при высоких давлениях в условиях контакта с аморфным инертным газом (аргоном и ксеноном). Поверхностное предплавление проявилось при моделировании для всех металлов. Для железа и молибдена температурный интервал T, в котором наблюдалось поверхностное предплавление, имел тот же порядок что и температурное расхождение между кривыми плавления полученными в статических и динамических экспериментах. Для алюминия поверхностное предплавление проявилось слабо Рис. 8: Слева: Сравнение экспериментальных и теоретических КП и области ПП для железа и алюминия. 1, 2, 3 кривые плавления измеренные в DAC-экспериментах; 4 и точки плавления по данным ударно-волновых экспериментов; 6 и 7 кривые плавления, вычисленные в (Laio et al., 2000) и (Alfe, 2009); 8 кривая плавления рассчитанная в данной работе и практически совпадающая с (Belonoshko et al., 2000). Заштрихованная область 12 зона предплавления для системы металл-аргон.

Справа: Расположение атомов (спроектированное на XY плоскость) в специальном моделировании деформации поверхности: синие атомы - аргон, красные атомы - железо.

Показана начальная геометрия (a) и геометрия после времени расчёта 30 нс (б).

по сравнению с погрешностями экспериментальных данных (рис. 8). Таким образом, поверхностное предплавление может быть одной из причин расхождения экспериментальных данных по определению кривой плавления железа и молибдена в области высоких давлений.

Для оценки свойства деформации поверхности металла при предплавлении был проведен расчет для системы железо-аргон с 285330 атомами (рис. 8). На протяжении времени расчета 30 нс наблюдалась эволюция первоначально сильно искривленной поверхности раздела компонент при условиях P = 167 ГПа и T = 4540 K (Tm-T = 120 K). При предплавлении поверхность действительно способна сильно деформироваться, устраняя первоначальные неровности, что должно приводить к изменению ее отражательной способности и, возможно, является одной из причин противоречий в экспериментальных данных полученных с использованием техники алмазных наковален (DAC, в данных экспериментах фиксируется преимущественно поверхностное плавление) по плавлению металлов при больших давлениях с данными ударно-волновых экспериментов.

Кавитация в жидкости и откольное разрушение. Жидкая фаза часто возникает в твердом теле в результате интенсивного импульсного энерговРис. 9: Слева: Сравнение кинетики роста с законом вязкостного роста. 1 - вязкостный рост (8); 2 - МД-моделирование; 3 - расчет по уравнению (7).

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

клада, как, например, в результате лазерной абляции или ударно-волновом воздействии. При этом возникает потребность изучить прочность жидкости при растяжении и предложить методы ее предсказания в достаточно широком диапазоне времен длительности импульса растяжения: от нескольких пикосекунд до микросекунд.

Время жизни метастабильной фазы в отдельном расчете определялось по резкому увеличению давления, который сопровождает появление и рост пузырька. В случае, когда появление критического зародыша в системе происходит случайным образом и среднее время жизни гомогенной фазы не зависит от момента начала наблюдения, распределение МД-траекторий по временам жизни гомогенной фазы будет подчиняться экспоненциальному закону N(t) = N0 exp (-t/ t ) (Скрипов и Коверда, 1984). Результаты МДмоделирования согласуются с данной моделью нуклеации.

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

Движение границы полости в несжимаемой жидкости описывается уравРис. 10: Слева: Зависимость откольной прочности гексана от скорости растяжения. 1 - расчет по модели NAG, 2 - расчет по нуклеационной модели, 3 - расчет по энергетическому критерию откола, 4 - экспериментальные данные (Уткин и др., 2003), 5 - прямое МД-моделирование, 6 - спинодаль.

Справа: Распределение полостей по размерам в момент откола t = tsp.

нением Рэлея-Плессета 3 4 2 P RR + 2 + + = -, (7) 2 R R где вязкость жидкости, поверхностное натяжение, плотность.

Из независимых МД-расчетов была определена зависимость вязкости от плотности для температуры T = 0.7 в том же диапазоне плотностей, в котором проводились МД-расчеты роста полостей.

В предположении, что рост полостей определяется, в первую очередь, вязкими свойствами среды, уравнение имеет очевидное решение (закон вязкостного роста) P R(t) = R0 exp - t. (8) 4 Сравнение кинетики роста полостей в МД-расчетах и соответствующих решений уравнения (7) показано на рисунке 9. Видно, что гидродинамический подход хорошо подходит для описания начального участка роста. На графике справа видно, что приближение вязкостного роста в случае ЛеннардДжонсовской жидкости неприменимо. Кроме того, выражение (8) дает неверную асимптотику решения. Точное решение уравнения (7) выходит асимп2P тотически на линейный рост R - при t . Вязкостный рост, 3 напротив, дает постоянное ускорение границы полости.

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

Затем система выводилась на равновесие при температуре T = 0.71, и начинался основной расчет. В ходе расчета объем МД-ячейки менялся с заданной скоростью. Расчеты проводились для различных скоростей растяжения. Для расчетов при высоких скоростях растяжения число атомов в ячейке равно 512 тыс., при низких 64 млн. На рис. 10 представлена полученная зависимость откольной прочности жидкости от скорости растяжения. На низких скоростях растяжения наблюдается слабая зависимость откольной прочности от скорости растяжения, на высоких выход откольной прочности на постоянное значение, связанный с достижением спинодали.

Было проанализировано распределение полостей по размерам в момент откола. Как видно из рис. 10, максимальную долю объема составляют полости большого радиуса. Они же дают основной вклад в скорость разрушения (скорость, с которой нарастает объем пустот). В то же время число этих полостей намного меньше числа маленьких полостей, зародившихся непосредственно перед отколом. Таким образом, разрушение происходит не за счет появления новых полостей непосредственно перед отколом, а за счет роста небольшого числа полостей, зародившихся несколько ранее.

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

На рис. 10 представлено сравнение расчетов откольной прочности по модели Унуклеация и ростФ (NAG). Для сравнения приведены расчеты откольной прочности по двум другим моделям. Также на рис. 10 нанесены результаты расчетов откольной прочности в прямом МД-моделировании и экспериментальные данные (Уткин и др., 2003). Хорошее совпадение результатов расчетов по модели Унуклеация и ростФ с экспериментальными данными и данными МД-экспериментов доказывает применимость предлагаемой модели к описанию разрушения.

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

Возмущения прямолинейности формы и равномерности движения дислокации в тепловом поле решетки вызывают индуцированное излучение дислокацией фононов. Движение дислокации в периодическом потенциале, свяРис. 11: Слева: Схема модели для расчета поведения дислокации в ГЦК кристалле под действием сдвига. Показаны только атомы, находящиеся в плоскости скольжения дислокаций. Темные атомы - дефект упаковки между двумя частичными дислокациями. Часть ячейки со свободными атомами M. Область с фиксированными атомами - F. Область атомов в виде абсолютно твердого тела R, которая двигается в направлении оси x с заданной скоростью vR.

Справа: Температурная зависимость коэффициента динамического торможения дислокаций B0(T ). Значения температуры и коэффициента торможения для алюминия из МД расчетов обезразмерены на температуру = 230K и величину B(). 1 - молекулярная динамика; 2 - аналитическая аппроксимация данных МД; 3 - эксперимент (Gorman et al., 1969); 4 - эксперимент (Hikata et al., 1970); 5 - теоретическая зависимость (Альшиц и Инденбом, 1975).

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

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

На рис. 11 представлена схема расчетной ячейки, соответствующая этой модели. Направления осей ориентированы соответственно одной из систем скольжения, типичной для ГЦК решетки, какой обладает алюминий, а именно [1 10](111). Вектор Бюргерса полной краевой дислокации, обеспечивающей такое скольжение, будет b = a/2[1 (направление совпадает с осью 10] x). Скольжение дислокации происходит в плоскости xz, т.е. (111). Линия краевой дислокации совпадает с осью z. Периодические граничные условия действуют вдоль осей x и z. Подвижные частицы расчетной ячейки составляют блок M (его размеры LbH L). Для создания сдвиговых напряжений xy используются несколько крайних атомных слоев, каждый состоящих из частиц, неподвижных друг относительно друга. Положение трех атомных слоев у нижней границы (F ) фиксируется на месте, в то время как слои у верхней границы (R) двигаются вдоль оси x под действием заданной внешней силы Fx, либо с постоянной заданной скоростью vx. Таким образом, в первом случае непосредственно контролируется сдвиговое напряжение, во втором деформация.

Для создания дислокации используется следующая процедура:

а) в идеальной кристаллической структуре удаляются две соседних атомных полуплоскости (1 10);

б) при небольшом сжатии кристалла происходит сближение атомных плоскостей, находящихся по разные стороны от места выреза, координаты атомов по осям y и z при этом фиксируются (т.е. движение атомов происходит вдоль оси x);

в) снятие ограничения на движение вдоль оси x для частиц блока M; УзамораживаниеФ частиц в блоках R и F, релаксация образовавшейся структуры методом минимизации потенциальной энергии;

г) релаксация с выводом на заданные температуру и давление.

Отметим важность третьего этапа (в) описанной процедуры, поскольку на нем происходит процесс расщепления полной дислокации, энергетически менее выгодной для ГЦК решетки, на две частичные краевые дислокации:

a a a [1 [2 + [1 10] 11] 21] 2 6 Методом МД были рассчитаны зависимости скорости дислокации v от сдвигового напряжения для нескольких температур, вплоть до температуры плавления. Четко выделяются два режима (особенно при низких температурах): линейный в области низких значений напряжения и режим асимптотического приближения скорости дислокации к поперечной скорости звука. Во всем исследованном диапазоне увеличение температуры приводит к уменьшению скорости дислокации, таким образом, движение дислокаций не требует термоактивации: имеет место динамический режим с фононным трением. Движение лимитируется перекачкой энергии от дислокации к элементарным возбуждениям в кристалле. Обычно при анализе зависимости v от ограничиваются линейным участком.

III [110] [111] Рис. 12: Слева: МД модель взаимодействия дислокации и препятствия.

Справа: Последовательные моменты времени процесса отрыва дислокации при температуре T = 300K (I) и T = 500K (II). Время между моментами - 4 пс (I) 40 пс (II). Показаны только атомы с дефектным локальным окружением.

Абсолютное значение коэффициента торможения дислокации при 300 К, определённое из обработки МД расчетов, составило приблизительно 0.9 10-5 Пас, что близко к экспериментальным значениям. По оценкам различных экспериментов оно варьируется от 0.5 10-5 Пас (затухание ультразвука, 1970 (Hikata et al., 1970)) и 2.6 10-5 Пас (подвижность дислокаций, 1969 (Gorman et al., 1969)) до 17 10-5 Пас (затухание ультразвука, 1966, см. обзор (Альшиц и Инденбом, 1975)). Значение коэффициента торможения из различных экспериментов систематически отличаются друг от друга, что может объясняться различием используемых методик и физической неэквивалентностью образцов. Поэтому вид температурной зависимости B(T ) различных экспериментов удобнее сравнивать именно в безразмерном виде (Альшиц и Инденбом, 1975), подбирая параметры , B() индивидуально (рис. 11).

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

Механизм взаимодействия дислокации в кристалле алюминия с включением меди показан на рис. 12). В ГЦК кристалле полная дислокация диссоциирует на две частичные, прохождение преципитата осуществляется в два этапа: отрыв ведущей дислокации, затем - ведомой. Но поскольку энергии дефекта упаковки, разделяющего две частичные дислокации в Al велика, то дислокациям не выгодно расходиться на большие расстояния и отрыв происГПа (iii) (ii) (i) Рис. 13: Слева: Зависимость пороговых напряжений, необходимых дислокации для прохождения препятствия, от радиуса препятствия r и расстояния между ними L. Кластер меди в Al. T = 300К (синие точки), T = 800К (красные точки). Сплошные линии - значения, рассчитанные по формуле (9) с A=2.9 Н/м и B=0.4 нм (синие), A=5 Н/м and B=0.нм (красные).

Справа: Температурная зависимость динамического предела текучести. 1 - экспериментальные результаты монокристаллов Al (Канель и др., 2007); 2 - экспериментальные результаты сплава Al D16T (Garkushin et al., 2008); 3 - оценка для монокристалла; 4 - критические напряжения для сплава Al с Cu (4%). Стрелки демонстрируют три эффекта: (i) - увеличение динамического предела текучести для монокристаллов с температурой, (ii) - упрочнение при введение включений, (iii) - разупрочнение сплава с температурой.

ходит практически одновременно.

На рис. 13 представлена зависимость cr от размера включения для нескольких расстояний между включениями и двух температур. На этом рисунке стоит отметить два факта: (i) увеличение размера включений приводит к увеличению критических напряжений, причем зависимость нелинейна; (ii) увеличение расстояния между включениями приводит к уменьшению значений критических напряжений. Это фактически связано с тем, что увеличение длины дислокации приводит к увеличению силы, действующей на всю ее длину.

Проведем сравнение полученных значений критических напряжений с предсказаниями теории. В рамках континуальной теории для такого механизма преодоления включения предложена следующая формула:

-A 2r 2r cr = log 1 +, (9) L B L где A = /2, L расстояние между включения, r их радиус, B свободный параметр. Для рассматриваемых размеров включений и расстояний между ними выполняется соотношение max(2r/L) 0.1 1. Поэтому формула может быть упрощена: cr = (A/L) [log(r) - log(B/2)].

Результаты МД расчетов пороговых напряжений прохождения дислокаций через препятствия представлены на рис. 13 для нескольких значений расстояния между препятствиями L и двух температур T. Для больших кластеров с ростом температуры наблюдается существенно меньшее по величине снижение предельного напряжения. Для кластеров малого размера, сравнимого с расстоянием между частичными дислокациями, и большим расстоянием между ними температурная зависимость оказывается обратной (см. точки 3-6 на рис. 13). Возрастание температуры сказывается также в увеличении времени задержки дислокаций на препятствиях, что является следствием роста величины фононного трения.

Величину напряжения течения при высокоскоростном нагружении можно определить, если считать, что скорость пластического деформирования определяется подвижностью дислокаций: = mbvD. Будем считать, что плотность подвижных дислокаций m постоянна в экспериментах, не зависит от температуры и условий деформирования. Тогда величина динамического предела текучести как функция температуры определяется температурной зависимостью коэффициента торможения дислокаций B(T ) в соответствии с формулой (Psakhie et al., 2007):

Y (T ) = B(T ) (10) mbПолученные таким образом значения можно сопоставить с данными о динамическом пределе текучести монокристаллического алюминия, полученными из ударно-волнового эксперимента (Канель и др., 2007), где скорости деформирования составляют приблизительно = 106 с-1 (см. точки 1 и кривую 3 на рис. 13). Точные значения плотности дислокаций в условиях ударно-волнового сжатия не известны, однако видно, что экспериментально полученная температурная зависимость хорошо описывается зависимостью динамического коэффициента торможения дислокаций от температуры B(T ), если величину плотности дислокаций принять равной m 7108 см-2.

Этот параметр является единственным свободным параметром и позволяет описать все точки экспериментальных данных ударно-волновых измерений. Различие на один - два порядка в сравнении с характерной плотностью подвижных дислокаций в монокристаллическом алюминии при нормальных условиях (106 - 107 см-2 по данным из (Gorman et al., 1969)) можно объяснить возрастанием плотности дислокаций при сжатии материала в ударной волне.

Пятая глава посвящена моделям откольного разрушения твердых тел.

Схема типичной экспериментальной постановки, позволяющей измерять Up Z [001] X [100] Y [010] Рис. 14: Схема типичного УВ эксперимента для определения динамической прочности материалов при импульсном деформировании.

прочность материала при отколе, представлена на рис. 14: ударник с некоторой скоростью налетает на мишень. В результате от плоскости соударения по ударнику и мишени начинают распространяться две ударные волны (УВ). За фронтом УВ активируется пластическое течение материала: происходит образование и движение дислокационных петель, дефектов упаковки.

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

Для исследования влияния свободной боковой поверхности была построена модель о.ц.к. Fe, в которой вдоль поперечного направления y не используются периодические граничные условия, благодаря чему возможна релаксация свободных поверхностей мишени в этом направлении (вдоль направления z использовались периодические граничные условия). На рис. 15 показан пример расчета для системы 2.7 млн. атомов в основной расчетной ячейке.

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

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

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

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

Для кристаллов с дислокационной подсистемой наблюдается существенное снижение откольной прочности в сравнении с бездефектными монокристал Рис. 15: Структурные превращения в модели ударник-мишень со свободными боковыми границами в направлении y (показана верхняя правая четверть расчетной ячейки, симметричная нижней части): а структура волны разрежения в момент начала образования несплошностей (атомы раскрашены соответственно координационному числу K (в координационной сфере радиуса rK = 3.1 ): светлые о.ц.к. структура (K = 14), темные плотноупакованная структура (K = 12), пунктиром показана плоскость симметрии модели); б увеличена область расчетной ячейки, близкая свободной боковой поверхности мишени (показаны только атомы с нарушенной симметрией локального окружения); в развитие откола через 5 пс после момента, соответствующего а и б. Скорость ударника uimp = 1400 м/с, полный размер системы 900x300х5 периодов решетки (2.7 млн. атомов).

ами. Для случая кристаллов с дефектами зависимость от скорости деформирования значительно более сильная и хорошо согласуется с зависимостью, полученной в ударно-волновых экспериментах (Kanel et al., 2001). Аналогичное поведение имеет место и для поликристаллов: зарождение полостей происходит на дефектной структуре межзеренных границ и зависимость откольной прочности от скорости деформирования является сильной.

МД расчеты позволяют детально проследить процесс роста пор на микроуровне, рассчитать кинетические характеристики их появления и скорости роста. Полученные таким образом данные о кинетике элементарных процессов составляют основу предлагаемой кинетической модели разрушения, относящейся к моделям Унуклеации и ростаФ несплошностей (Barbee et al., 1972). В результате такая кинетическая модель позволяет перейти от атомистических механизмов разрушения к описанию процессов образования и развития несплошностей в материале на макроскопическом уровне.

На рис. 17 представлен график зависимости рассчитанной величины откольной прочности psp в динамическом режиме деформирования (растяжение с постоянной скоростью) от скорости растяжения . В полулогариф мическом масштабе зависимость приблизительно линейная, поскольку разрушение определяется процессом гомогенной нуклеации полостей, скорость которого экспоненциально зависит от степени деформации.

В рассмотренном на правой части рис. 17 случае высоких температур реРис. 16: Эволюция микроструктуры монокристалла в процессах: предварительного одноосного сжатия (a), релаксации к нулевым напряжениям (b) и гидростатического растяжения на (c,d). Использован потенциал (Mishin et al., 1999). Число атомов в системе 2 106. Атомы с локальным окружением, соответствующим г.ц.к. решетке, не показаны.

Вакансии отмечены на рис. (b,c) кругами, а дислокации и дефекты упаковки - прямоугольниками.

зультаты расчета находятся в хорошем согласии как с экспериментальными данными при более низких скоростях деформирования, так и с прямым МД моделированием в области высоких скоростей. Прямое МД моделирование представляет собой всестороннее растяжение монокристаллического Al при той же температуре. Изменение откольной прочности при изменении скорости деформирования тесно связано с зависимостью частоты зарождения полостей от напряжения, которая меняется, если имеет место плавление. Таким образом, при высоких температурах плавление важно для оценки откольной sp, GPa sp, GPa - - 8 - - NAG - , s-, s-105 107 109 10105 107 109 10Рис. 17: Слева: Зависимость откольной прочности от скорости деформирования. Зависимость из МД расчетов гидростатического растяжения при постоянной температуре T = 300 K: 1 - бездефектный кристалл, потенциал (Mishin et al., 1999). 2,3 - модель с дефектами, потенциалы (Liu et al., 2004; Mishin et al., 1999) соответственно; 4 - экспериментальные данные (Kanel et al., 2001). 5 - аналитическая аппроксимация.

Справа: Сравнение данных по откольной прочности монокристаллического Al при T 900 K: прямое МД моделирование (данная работа); кинетическая модель NAG (данная работа); экспериментальные данные (Kanel et al., 2004).

прочности: оно приводит к более сильной зависимости откольной прочности от скорости растяжения. Хорошая точность описания разрушения вблизи температуры плавления в рамках модели Унуклеация и ростФ связана с близостью к жидкой фазе, где эта модель работает надежно (см. главу 3).

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

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

На рис. 18 показаны три фрагмента со структурой атомов, соответствующие трем различным стадиям эволюции облученного металла: первоначальный трек, первоначальное распределение дефектов сразу после кристаллиРис. 18: Три фрагмента эволюции структуры о.ц.к. Mo после пролета и остановки иона Xe. Температура системы 500 K. Верхний рисунок направление первоначального трека атома ксенона (отмечено красной стрелкой). Средний рисунок после времени эволюции 10 пс. Нижний рисунок после времени эволюции 4 нс. Показаны только атомы с потенциальной энергией выше -6.5 эВ. На вставках приведено расположение всех атомов спроецированное на XY-плоскость: а виден один подвижный дефект из четырех SIA и один неподвижный кластер из пяти SIA; 2 подвижных di-SIA (по разному ориентированных) и один подвижный кластер из четырех SIA.

зации трека (10 пс после облучения) и конечное распределение дефектов после времени эволюции 4 нс. На данном рисунке (за исключением вставок) показаны только атомы с потенциальной энергией выше -6.5 эВ. Т.к. дефект межузельного атома создает довольно сильною деформацию вдоль оси <111>, то мы видим только сами межузельные атомы и несколько ближайших соседей вдоль оси <111>. Хорошо видно, что в первоначальный момент времени присутствуют только единичные дефекты. Цепочка из 5 атомов на рисунке представляет собой одиночный межузельный атом (SIA). Вакансии на рисунке не показаны, но можно сказать, что они практически не эволюционируют в течении времени расчета и сохраняют первоначальное положение (их расположение совпадает с треком ксенона). В течение первых пс после облучения выкристаллизовываются около 100-150 SIA и вакансий. Однако большинство дефектов практически сразу рекомбинируют с вакансиями.

После времени эволюции системы равному 4 нс формируются и выживают только достаточно большие кластеры из дефектов.

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

Исходная структура МД для расчета скорости кластеризации SIA в двойные кластеры (di-SIA) и рекомбинации соответствует о.ц.к. кристаллу молибдена с добавленными точечными дефектами. В течении эволюции дефектов во время МД расчета новых дефектов не создается. Первоначально вводится два типа дефектов: межузельные атомы (SIA) в форме краудионов и вакансии (в случае расчетов рекомбинации). Начальное пространственное распределение дефектов является однородным для дефектов обоих типов.

Краудионы имеют четыре эквивалентных направления в элементарной ячейке: [111], [ [1 и [11 Вероятность задания краудиона каждой ори111], 11] 1].

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

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

Изменение со временем концентрации межузельных атомов ci и вакансий cv описывается следующей системой кинетических уравнений:

dci dcv 2 2 2 2 2 = G-kvDici-ki Dvcv-kisDici, = G-kvDici-ki Dvcv-kvsDvcv, (11) dt dt где G - скорость производства дефектов в cна/c, Di, Dv и ki, kv - коэффициенты диффузии и скорости стока межузельных атомов и вакансий соответственно. Второе и третье слагаемые в правой части (11) описывают рекомбинацию в результате диффузии межузельных атомов к вакансиям и вакансий 2 к межузельным атомам (поэтому ki и kv это силы вакансий, выступающих в качестве стока для межузельных атомов, и наоборот). Последние слагаемые соответствуют кластеризации межузельных атомов и вакансий.

Вообще говоря, подобная системы уравнений может решаться для конкретных задач путем численного интегрирования или, в некоторых случаях, аналитически. Наша МД модель соответствует случаю G = 0. Коэффициент диффузии вакансий Dv намного меньше, чем коэффициент диффузии междуузлий Di во всем рассмотренном температурном диапазоне. Для упрощения анализа мы будем считать диффузию трехмерной. Будет учитываться зависимость коэффициентов от температуры. В итоге, система (11) приобретает более простой вид:

dci dcv = -rcicv - iicici, = -rcicv, (12) dt dt где r и ii имеют смысл скоростей бинарных реакций.

Рис. 19 иллюстрирует изменение со временем числа межузельных атомов и кластеров di-SIA в системе, не содержащей вакансий. В течении первых пс система приходит в квазиравновесное состояние в постоянным температурой и давлением. Близко расположенные межузельные атомы взаимодейсвуют в процессе атомной динамики и образуют кластеры. Метод диагностики дефектов и их подсчета проверяется сравнением полного числа дефектов Ntotal = NSIA +Ndi-SIA и начального числа дефектов N0. Отличие этих величин показывает число образовавшихся тройных и больших кластеров и, тем самым, показывает, с какой точность выполняется предположение о том, что число кластеров, больших di-SIA, пренебрежимо мало. В главе приводятся результаты аналогичных расчетов для рекомбинации SIA и вакансий.

Таким образом, предложена концепция многомасштабной модели для расчета характеристик эволюции дефектов на примере в молибдена. Данная концепция сочетает построение межатомного потенциала на основе ab initio расчетов, моделирования эволюции ансамблей дефектов на уровне классической МД, подбор начальных условий МД моделей для расчета конкретных enabling of our MDEq. withcan kinetic theory resultsthen the transition fromcomparisondiffusion. (1) be simplified:to be 3D to 1D made in the coordinates 1/ci 1/ci0 vs. t, in order to obtain the rate dci constant.

arcicv aiicici;

dt 2 Fig. 2 shows the evolution of defect concentrations obtained by dcv atathree:initial concentrations of SIAs at 300 K. The dependencicv r dt MD cies are well approximated by the linear equation with the slope Here ar which are the rateto be nearly terms of theof the initial cona constants inindependent binary reacpffiffiffiii aii,and turned out 2=tions; aii 2Di The resulting rates are for association of3two SIAs centration.=X is the rate constantpresented in Fig. for different into one di-SIA. Eq. (2) can be analytically in order to initial concentrations andsolved temperatures. One cancomfor two conclude pare with MD results.

that the reaction is second order at both temperatures, although Fig. 3. Dependence of rate constant aii on the initial concentration for two temperatures T: 1Ц300 K, 2Ц1000 K. The system size is 2 million atoms.

the motion is strongly 1D at T = 300 K and mostly 3D at T = 1000 K.

Fig. 4a shows the time dependence of the concentrations of mo3.3. MD calculations of kinetic coefficients bile and sessile clusters during clustering of SIAs. An interesting feature of this plot is that the number of sessileof saturates clusters Fig. 1 illustrates the evolution of the number interstitial ably transform to the mobile configuration and reach a quasi-equiabout whereas the number of mobile clusters is still atomsafter di-SIAs30 a without vacancies. During the first and inps,system Fig. 2. Dependence of the SIA number on time for different initial concentrations librium concentration.

growing. This situation means that the sessilewhere most prob- clusters 10 ps the system reaches quasi-equilibrium state, the tem- ci0: 1Ц1.5 10 4,, 3Ц3 10 4. T = 300 K.

Fig. 4b2Ц2 10 that the relative number of sessile cluster beshows perature and pressure are equilibrated and the very close intersticomes less at T = 1000 K, which characterizes thermally activated tials interact by forming clusters. The diagnostics is verified by the transformation of the sessile clusters into mobile ones. Based on total number of recognized defects, Ntotal, equal to NSIA +Ndi-SIA. If the temperature dependence of the quasi-equilibrium concentraNtotal equals its initial value N0, the diagnostics accurately recogtion of sessile clusters, one can estimate activation energy of trannizes all defect structures. After a long calculation time the Ntotal sition into mobile clusters as 0.2 eV.

changes, indicating that three SIAs and larger clusters are formed.

Integration of Eq. (2) with cv = 0 gives the following solution:

3.4. Determination of rate constants 1=ci 1=ci0at; 3 ii The recombination rates of SIA and VAC were determined from MD simulation of the systems, where the concentration of vacanenabling comparison of our MD results with the kinetic theory to be cies is much larger than the concentration of interstitials. It premade in the coordinates 1/ci 1/ci0 vs. t, in order to obtain the rate sents the possibility of neglecting the change of vacancies constant.

concentration: cv(t) cv0. In this case Eq. (2) results in the followFig. 2 shows the evolution of defect concentrations obtained by ing simple solution:

MD at three initial concentrations of SIAs at 300 K. The dependencies are well approximated by the linear equation with the slope cici0 exp aivcv0t; 4 aii, which turned out to be nearly independent of the initial concentration. The resulting rates are presented in Fig. 3 for different The concentrations of interstitials were fitted according to Eq. (4) at initial concentrations and for two temperatures. One can conclude five temperatures. Fitting procedures were conducted at several inithat the reaction isof the number of SIAs and di-SIAs, and total number of recognized Fig. 3. Dependence of rate of vacancies, ensure that the results do second order at both temperatures, although tial concentrations constant aii ontothe initial concentration for twonot Fig. 1. Evolution temperatures 1Ц300 2Ц1000 K. The system is 2 million Рис. 19: Слева: Эволюция числа SIA, двойных кластеров (di-SIA) и полное количество the motion is stronglysize at T = 300 K and mostly 3D at T = 1000 K.

1Dis defects. The system 2 million atoms, and T = 1000 K. dependT:on theK, concentrationsize vacanciesatoms. 5). One can initial of (Fig.

Fig. 4a shows the time dependence of the concentrations of moраспознанных дефектов. Размерinteresting 2 миллиона атомов, T = 1000 K.

системы bile and sessile clusters during clustering of SIAs. An Please cite this article in press as: Z. Insepov et al., J. Nucl. Mater. (2011), doi:10.1016/j.jnucmat.2011.08.0Справа: Зависимость константы скорости feature of this plot is that the number of sessile clusters saturates ii от начальной концентрации для двух темably transform to the mobile configuration and reach a quasi-equiafter about 30 ps, the number mobile clusters is still ператур:whereas К, 2 Цof K. Размер системы 2 миллиона атомов.

1 - 300 10librium concentration.

growing. This situation means that the sessile clusters most probFig. 4b shows that the relative number of sessile cluster becomes less at T = 1000 K, which characterizes thermally activated transformation of the sessile clusters into mobile ones. Based on кинетических констант, методы проверки независимости результатов расчеthe temperature dependence of the quasi-equilibrium concentraтов от выбора параметров моделей, не связанных с заданием макроскопичеtion of sessile clusters, one can estimate activation energy of transition into mobile clusters as 0.2 eV.

ских параметров.

3.4.

В седьмой главе описаны моделиDetermination of rate constantsразогретого плотнонеравновесного The recombination rates of SIA and VAC were determined from го состояния (warm dense matter) при абляции диэлектриков и металлов на MD simulation of the systems, where the concentration of vacanпримерах фторида лития и золота.

cies is much larger than the concentration of interstitials. It presents the possibility Современные технологии дают возможность(t) c of neglecting the быстрых энербеспрецедентноchange of vacancies concentration: cv v0. In this case Eq. (2) results in the followговкладов в конденсированную фазу,ing simple solution: к новым типам задач для что приводит cici0 exp aivcv0t; 4 теории и компьютерного моделирования. Электромагнитное поле создает The concentrations of interstitials сильные возбуждения в электронной системе твердогоwere fitted according to Eq. (4) at к тела и приводит five temperatures. Fitting procedures were conducted at several initial concentrations of vacancies, ensure that the results do not Fig. 1. Evolution of the number of SIAs and di-SIAs, and total number of recognized двухтемпературному разогретому плотному состоянию.to Поскольку5). One can соответdefects. The system size is 2 million atoms, and T = 1000 K. depend on the initial concentration of vacancies (Fig.

ствующее изменение в электронном распределении может быть практически Please cite this article in press as: Z. Insepov et al., J. Nucl. Mater. (2011), doi:10.1016/j.jnucmat.2011.08.0мгновенным в сравнении с характерными временами движения ядер, происходит быстрое изменение эффективного ландшафта поверхности потенциальной энергии кристаллической решетки и сил действующих н атомы. Это может приводить к потере устойчивости кристаллической решетке и ее аморфизации (т.н. нетермическому плавлению) до того, как энергия переходит от электронов в решетке в результате электрон-фононного взаимодействия.

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

Исходная теория функционала электронной плотности, направленная на расчет основанного состояния многоэлектронной системы во внешнем потенциале, была обобщена на случай конечных температур путем рассмотрения электронных возмущений в усредненном статистическом смысле (Mermin, 1965).

Этот подход был использован в работе (Alavi et al., 1994) и последующих -20 -10 0 10 20 30 E, eV ne(r)=-0.F ne(r)=+0.Li X L Рис. 20: Слева: Элементарная ячейка г.ц.к. кристалла LiF. Показаны поверхности постоянного значения разницы объемных распределений электронной плотности при различных электронных температурах ne = ne(r)|3.2eV - ne(r)|0.0eV : ne = -0.04 (сферическая поверхность вокруг атома) и ne = +0.04 (две звездообразных области между Li и F).

Справа: Акустические ветви дисперсионных зависимостей фононов в г.ц.к. решетке LiF при различных электронных температурах: Te = 2.0 эВ (сплошные линии) и Te = 3.2 эВ (пунктирные линии).

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

Фторид лития. Эксперименты на фемтосекундном рентгеновском лазере по сверхбыстрому энерговкладу в LiF (Faenov et al., 2009) показали аномально низкие пороги абляции, свидетельствующие о возможном сильном изменении механических свойств фторида лития в состоянии разогретого плотного вещества. Цель исследования, представленного в этой главе, заключается в использовании первопринципных квантово-механических расчетов для описания изменений в межатомных взаимодействиях в кристалле LiF в данном переходном неравновесном состоянии.

Расчеты были выполнены с помощью пакета VASP. 2s электроны атомов лития и 2s и 2p электроны атомов фтора рассматриваются в качестве валентных. Используются PAW псевдопотенциалы для описания совместного потенциала ядер и электронов кора (Kresse and Furthmuller, 1996). Обменнокорреляционная компонента полной энергии системы описывается в градиентном приближении (GGA). Для расчета акустических ветвей дисперсионных зависимостей фононов используется метод малых смещений (Alf, 2009).

Увеличение электронной температуры приводит к изменению плотности электронных состояний в кристалле. Энергии s-зон понижаются в сравнении DOS, states/eV Frequency, THz с нормальными условиями. Соответствующее перераспределение электронной плотности в элементарной ячейке показано на рис. 20. Видно, что электронная плотность s-электронов, локализованная при Te = 0 вокруг атома фтора, частично переходит в две звездоподобных области между атомами Li и F. Это свидетельствует о том, что имеющая место при нормальных условиях сильная ионная связь LiF приобретает ковалентный характер в результате возбуждения электронной подсистемы.

Сопутствующее изменение эффективного межионного потенциала должно приводить к изменениям в механической устойчивости г.ц.к. решетки. Для количественного описания этих изменений были проведены расчеты дисперсии акустических колебаний в решетке. Результаты для Te = 0 находятся в хорошем согласии с недевними подробными расчетами (Evarestov and Losev, 2009). Рис. 20 показывает, что решетка остается устойчивой при электронной температуре 2.0 эВ, однако при температуре 3.2 эВ появляется мягкая акустическая мода, что говорит о потере механической устойчивости исходной г.ц.к. структуры и возможности нетермического плавления. Описание г.ц.к. LiF в градиентном (GGA) приближении занижает величину запрещенной зоны широкозонного диэлектрика LiF. Поэтому приведенные значения электронной температуры имеют более качественный, чем строго количественный характер.

Возбуждение электронной подсистемы до температур более Te 3 эВ приводит к потере механической устойчивости г.ц.к. кристалла LiF, что проявляется в виде возникновения мягкой моды в спектре акустических фононов и, по-видимому, должно приводить к нетермическому плавлению. Соответствующее перераспределение электронной плотности подразумевает, что исходная сильная ионная связь LiF приобретает более ковалентный характер с ростом электронной температуры.

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

В то же время, механизм лазерной абляции металлов остается не вполне ясным. На это указывают как противоречия в экспериментальных данных, так и отсутствие хорошего согласия между экспериментами и предсказаниями теории. Одним из наиболее изучаемых металлов является золото, которое характеризуется большим, чем остальные металлы, временем релаксации электронной температуры Te к температуре ионов Ti (Ernstorfer et al., 2009;

Lin et al., 2008). Поэтому золото является удобным металлом для исследования эффекта двухтемпературности (Te Ti) и построения общей модели лазерной абляции металлов. Важной непосредственно измеряемой характеристикой абляции является глубина кратера d, другой важной характеристиРис. 21: Потенциальные функции EAM-потенциала для золота для трех различных температур: синие сплошные линии при Te = 0.1 эВ; коричневые штриховые линии при Te = 3 эВ; красные штрих-пунктирные линии при Te = 6 эВ.

кой является минимальная пороговая величина поглощенного флюенса Fabs, необходимая для начала абляции.

Атомистическая модель двухтемпературной системы, предложенная в работах (Ivanov and Zhigilei, 2003; Schafer and Urbassek, 2002), заключалась в том, что электронная подсистема моделировалась сплошной средой. Эволюция ионной подсистемы описывалась с помощью метода молекулярной динамики. Передача энергии от электронной к ионной подсистеме выполнялась путем введения термостата в МД модель. Такой подход позволяет объяснить многие особенности лазерной абляции (Demaske et al., 2010; Ivanov and Zhigilei, 2003), но дает завышение порогового порогового поглощенного флюенса, по крайней мере, для Cu (Hu et al., 2010; Schafer and Urbassek, 2002) и для Au (Demaske et al., 2010). В этих исследованиях абляция связана с формированием и распространением ударной волны и волны разгрузки вглубь металла. Указанное противоречие указывает на сложности описания абляции при энерговкладах вблизи порогового флюенса.

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

Для создания ETD-потенциала по результатам ab initio расчетов использовалась методика, аналогичная описанной в главе 2. На рис. 21 представлены найденные потенциальные функции (r), (r) и () для разных значений Te. Потенциалы были сконструированы таким образом, чтобы функция Рис. 22: Слева: Зависимость давления золота от Te: 1 полное давление вычисленное с помощью ab initio расчетов (использовался VASP код); 2 давление в системе, описанной ETD-потенциалом, восстановленное по теореме вириала в молекулярно-динамическом deloc расчете. Разница между зависимостями определяет Pe.

Справа: Изотермы золота при различных Ti и Te (полученные с использованием ETDdeloc потенциала и без учета Pe ): 1 Ti = Te = 0.01 эВ; 2 Ti = Te = 0.2 эВ (жидкость); 3 Ti = 0.01 эВ, Te = 3 эВ. Стрелкой схематично показана эволюция состояния вблизи поверхности при УкороткойФ абляции в ходе расчета с = 0.1 пс, l = 6 нм и deloc Fabs = 121 мДж/см2 (с учетом Pe ): а состояние, соответствующее нормальным условиям; b быстрое нагревание электронной подсистемы лазерным импульсом; c релаксация вследствие расширения системы и уменьшения Te; d конечное состояние, соответствующее P - 5.6 ГПа, в котором происходит разрушение золота.

() была одинакова для всех Te. Все три потенциала воспроизводят силы, рассчитанные ab initio, с точностью 15-20 процентов, что является хорошей точностью для данного метода (Brommer and Gahler, 2007).

Следующий этап разработки заключался в создании параметрической зависимости функций (r), (r) от Te. Это было сделано путем представления указанных функций в виде квадратичного полинома от Te:

(r, Te) = 0(r) + 1(r)Te + 2(r)Te, (13) (r, Te) = 0(r) + 1(r)Te + 2(r)Te. (14) Данные функциональные зависимости вместе с функцией () представляют собой разработанный ETD-потенциал, который, в данном случае, определяется семью независимыми функциями.

Разработанный ETD-потенциал при низком значении Te хорошо воспроизводит термодинамические и механические свойства золота, известные из эксперимента (уравнение состояния, кривую плавления, упругие константы и др.) Хорошее согласие результатов моделирования с экспериментальными данными демонстрирует надежность использования ETD-потенциала не только для исследования двухтемпературного состояния, но и для описания равновесных свойств золота, когда Te = Ti.

Одним из эффектов, который воспроизводит ETD-потенциал является увеличение температуры плавления Tmelt, i при нагреве электронной подсистемы (Ernstorfer et al., 2009; Recoules et al., 2006). Температура плавления рассчитывалась с помощью метода двухфазного моделирования.

На рис. 22 показаны зависимости полного давления P (вычисленного в loc DFT расчете) и давления Pe (вычисленного с помощью МД с ETD-потенциалом) от значения Te. Разница между указанными зависимостями опредеdeloc ляется величиной Pe (Te), которая объясняется вкладом УделокализованныхФ электронов. Таким образом, другим эффектом, который воспроизводит ETD-потенциал, является увеличение давления УлокализованныхФ электроloc нов Pe при росте Te. При увеличении Te с 0 до 6 эВ давление в системе, описанной ETD-потенциалом, возрастает с 0 до 100 ГПа. Вместе с давdeloc лением УделокализованныхФ электронов Pe суммарное давление в золоте достигает 200 ГПа при Te = 6 эВ.

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

Для моделирования процесса лазерной абляции золота использовалась атомистическая двухтемпературная модель с ETD-потенциалом. Модель была реализована как развитие МД кода LAMMPS (Plimpton, 1995). Электронная подсистема описывалась как непрерывная среда с температурой Te.

Эволюция Te определялась выражением:

Te I(t)e-x/l Ce = (keTe) - G(Te - Ti) -, (15) t l где Ce электронная теплоемкость (все основные расчеты проводились с использованием зависимости Ce(Te) из работы (Lin et al., 2008), однако для верификации стабильности модели было проведено несколько расчетов с зависимости Ce(Te) из работы (Chen et al., 2006)), электронная теплопроводность ke и константа электрон-ионного взаимодействия G были взяты из работы (Ivanov and Zhigilei, 2003), I эффективная интенсивность лазерного импульса с длительностью , l глубина поглощения. В работе использовалось два набора и l: для моделирования оптических фемтосекундных импульсов = 0.1 пс и l = 6 нм; для моделирования рентгеновских пикосекундных импульсов = 7 пс и l = 18 нм.

Ионная подсистема описывалась как молекулярно-динамическая модель с Рис. 23: Фрагменты моделирования с расположением атомов золота в различных расчетах: (a) и (b) при = 7 ps с Fabs = 65 мДж/см2 и Fabs = 130 мДж/см2; (c) и (d) при = 0.1 ps с Fabs = 65 мДж/см2 и Fabs = 152 мДж/см2. Показана проекция атомов на xy-плоскость в расчетной ячейке.

уравнением движения:

deloc dvi (Pe ) mi = Fi(Te) + FLang(Te, Ti) -, (16) i dt ni deloc Рис. 24: Слева: Профили давления (без учета Pe ) вдоль x в разные моменты времени в ходе расчета с = 0.1 пс, l = 6 нм и Fabs = 152 мДж/см2: 1 t = 0.1 пс; 2 t = 3.5 пс; t = 20 пс (максимум давления при отрицательном значении x соответствует отколовшейся части металла); 4 t = 70 пс.

Справа: Зависимость d(Fabs) при лазерной абляции золота для = 0.1 пс: 1 и 2 значения полученные экспериментально с оптическими субпикосекундными импульсами (Inogamov et al., 2008; Vorobyev and Guo, 2005); 3 результат атомистического расчета без учета электронного давления (Demaske et al., 2010); 4 результаты моделирования в данной работе (с зависимостью Ce(Te) из работы (Lin et al., 2008); 5 то же, что и 4, но с зависимостью Ce(Te) из (Chen et al., 2006).

где сила действующая на ион складывалась из трех составляющих. Первое слагаемое является силой ETD-потенциала, зависящей от локального значения Te. Второе слагаемое является случайной силой Ланжевеновского термостата, обеспечивающего передачу тепловой энергии от электронной подсистемы к ионной. Последнее слагаемое представляет собой силу увлечения идеального газа делокализованных электронов (Chen et al., 2006; Gan and Chen, 2009), где ni ионная концентрация.

В данной работе показано, что УкороткаяФ абляция приводит к формированию кратеров на поверхности металла с глубиной порядка десяти нанометров и требует меньшего энерговклада, чем другие типы абляции (рис. 24).

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

механического расширения и уменьшения Te. Гипотеза об этом механизме позволяет согласовать экспериментальные результаты по абляции золота.

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

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

Классический масштаб:

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

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

3. Созданы модели динамики одиночной краевой дислокации в алюминии и ее взаимодействия с нановключениями меди. Рассчитаны температурная зависимость коэффициента подвижности дислокации и критические напряжения преодоления нановключений. Пластичность металлов при высокоскоростной деформации объяснена на основе результатов атомистических моделей.

Кинетический масштаб:

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

5. Определен столкновительно-лимитированный характер кинетики движения фронта плавления в металлах. Выявлен эффект предплавления металлов в контакте с внешней средой при высоких давлениях.

6. С использованием разработанного наиболее точного на данный момент межатомного потенциала для системы Mo-Xe выполнено моделирование трека иона Xe в о.ц.к. Mo и каскадных повреждений. Предложен метод расчета констант реакций взаимодействия различных типов радиационных дефектов в металлах. Проведены расчеты скоростей кластеризации одиночных межузельных атомов и их рекомбинации с вакансиями.

Квантовый масштаб:

7. Развит метод создания межатомных потенциалов в рамках модели потенциала погруженного атома для металлов на основе первопринципных расчетов. Создан межатомный потенциал для молибдена, превосходящий аналоги по точности описания его механических и термодинамических свойств, а также энергетической иерархии точечных дефектов.

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

8. На основе первопринципных расчетов дисперсии фононов показана потеря устойчивости кристаллической решетки фторида лития в разогретом плотном состоянии.

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

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

Т. 119, № 5. С. 1011Ц1020.

[2] Norman G. E., 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] Норман Г. Э., Стегайлов В. В. Гомогенная нуклеация в перегретом кристалле.

Молекулярно-динамический расчет // Доклады Академии Наук. 2002. Т. 386, № 3.

С. 328Ц332. (представлено академиком В.П.Скриповым) [4] Norman G. E., Stegailov V. V. Simulation of Ideal Crystal Superheating and Decay // Molecular Simulation 2004. Vol. 30, no. 9. Pp. 397Ц406.

[5] Kuksin A. Y., Morozov I. V., Norman G. E., Stegailov V. V., Valuev I. A. Standards for molecular dynamics modelling and simulation of relaxation // Molecular Simulaton. 2005. Vol. 31, no. 14-15. Pp. 1005Ц1017.

[6] Stegailov V. Homogeneous and heterogeneous mechanisms of superheated solid melting and decay // Computer Physics Communications. 2005. Vol. 169, no. 1-3. Pp. 247Ц250.

[7] Норман Г. Э., Стегайлов В. В., Янилкин А. В. Разрушение кристаллического железа при высокоскоростном растяжении. Молекулярно-динамический расчет // Доклады Академии Наук. 2005. Т. 404, № 6. С. 757Ц761. (представлено академиком В.Е.Фортовым) [8] Kuksin A., Norman G., Stegailov V., Yanilkin A. Surface melting of superheated crystals. Atomistic simulation study // Computer Physics Communications. 2007. Vol. 177, no. 1-2.

Pp. 34Ц37.

[9] Insepov Z., Hassanein A., Bazhirov T. T., Norman G. E., Stegailov V. V. Molecular Dynamics Simulations of Bubble Formation and Cavitation in Liquid Metals // Fusion Science and Technology. 2007. Vol. 52. Pp. 885Ц889.

[10] Стегайлов В. В., Янилкин А. В. Структурные превращения в монокристаллическом железе при ударно-волновом сжатии и растяжении. Исследование методом молекулярной динамики // ЖЭТФ. 2007. Т. 131, № 6. С. 1064Ц1072.

[11] Куксин А. Ю., Норман Г. Э., Стегайлов В. В. Фазовая диаграмма и спинодальный распад метастабильных состояний Леннард-Джонсовской системы // Теплофизика Высоких Температур. 2007. Т. 45, № 1. С. 43Ц55.

[12] Норман Г. Э., Стегайлов В. В., Янилкин А. В. Моделирование высокоскоростного растяжения кристаллического железа методом молекулярной динамики // Теплофизика Высоких Температур. 2007. Т. 45, № 2. С. 193Ц202.

[13] Куксин А. Ю., Стегайлов В., Янилкин А. В. Молекулярно-динамическое моделирование динамики краевой дислокации в алюминии // Доклады Академии Наук. 2008. Т. 420, № 4. С. 467Ц471. (представлено академиком Ю.А.Осипьяном) [14] Куксин А. Ю., Стегайлов В. В., Янилкин А. В. Атомистическое моделирование пластичноcти и разрушения нанокристаллической меди при высокоскоростном растяжении // Физика Твердого Тела. 2008. Т. 50, № 11. С. 1984Ц1990.

[15] Bazhirov T. T., Norman G. E., Stegailov V. V. Cavitation in liquid metals under negative pressures. Molecular dynamics modeling and simulation. // Journal of Physics: Condensed matter.

2008. Vol. 20, no. 11. P. 114113.

[16] Стариков С. В., Стегайлов В. В. Молекулярно-динамическое моделирование предплавления железа при высоком давлении // Доклады Академии Наук. 2009. Т. 424, № 1.

С. 31Ц35. (представлено академиком В.Е.Фортовым) [17] Kuksin A., Norman G., Stegailov V., Yanilkin A., Zhilyaev P. Dynamic fracture kinetics, influence of temperature and microstructure in the atomistic model of aluminum // International Journal of Fracture. 2009. Vol. 162, no. 1-2. Pp. 127Ц136.

[18] Starikov S. V., Stegailov V. V. Atomistic simulation of the premelting of iron and aluminum:

Implications for high-pressure melting-curve measurements // Physical Review B. 2009.

Vol. 80, no. 22. Pp. 20Ц23.

[19] Stegailov V. Stability of LiF Crystal in the Warm Dense Matter State // Contributions to Plasma Physics. 2010. Vol. 50, no. 1. Pp. 31Ц34.

[20] Жиляев П. А., Куксин А. Ю., Стегайлов В. В., Янилкин А. В. Влияние пластической деформации на разрушение монокристалла алюминия при ударно-волновом нагружении // Физика Твердого Тела. 2010. Т. 52, № 8. С. 1508Ц1512.

[21] Куксин А. Ю., Норман Г. Э., Писарев В. В., Стегайлов В. В., Янилкин А. В. Кинетическая модель разрушения простых жидкостей // Теплофизика Высоких Температур. 2010.

Т. 48, № 4. С. 536Ц543.

[22] Kuksin A., Norman G., Pisarev V., Stegailov V., Yanilkin A. Theory and molecular dynamics modeling of spall fracture in liquids // Physical Review B. 2010. Vol. 82, no. 17. Pp. 1Ц10.

[23] Янилкин А. В., Жиляев П. А., Куксин А. Ю., Норман Г. Э., Писарев В. В., Стегайлов В. В.

Применение суперкомпьютеров для молекулярно-динамического моделирования процессов в конденсированных средах // Вычислительные Методы и Программирование. 2010.

Т. 11. С. 111Ц116.

[24] Norman G. E., Skobelev I. Y., Stegailov V. V. Excited States of Warm Dense Matter // Contributions to Plasma Physics. 2011. Vol. 51, no. 5. Pp. 411Ц418.

[25] Starikov S., Insepov Z., Rest J., Kuksin A. Y., Norman G., Stegailov V., Yanilkin A. Radiationinduced damage and evolution of defects in Mo // Physical Review B. 2011. Vol. 84, no. 10.

Pp. 1Ц8.

[26] Derivation of kinetic coefficients by atomistic methods for studying defect behavior in Mo / Z. Insepov, J. Rest, A. Yacout, A. Kuksin, G. Norman, V. Stegailov, S. Starikov, A. Yanilkin // Journal of Nuclear Materials. 2011. DOI: 10.1016/j.jnucmat.2011.08.0[27] Лазерная абляция золота: эксперимент и атомистическое моделирование / С. Стариков, В. Стегайлов, Г. Норман, В. Фортов, М. Ишино, М. Танака, Н. Хасегава и др. // Письма в ЖЭТФ. 2011. Т. 93, № 11. С. 719Ц725.

итература В. И. Альшиц и В. Л. Инденбом. Динамическое торможение дислокаций. Успехи физических наук, 115:3, 1975.

Г. И. Канель, В. Е. Фортов, и С. В. Разоренов. Ударные волны в физике конденсированного состояния. Успехи физических наук, 177:809, 2007.

В. П. Скрипов и В. П. Коверда. Спонтанная кристаллизация переохлажденных жидкостей.

Наука, М., 1984.

А. В. Уткин, В. А. Сосиков, и А. А. Богач. Импульсное растяжение гексана и глицерина при ударно-волновом воздействии. ПМТФ, 44(2):27Ц33, 2003.

A. Alavi, J. Kohanoff, M. Parrinello, and D. Frenkel. Ab initio molecular dynamics with exited electrons. Phys. Rev. Lett., 73:2599Ц2602, 1994.

D. Alfe. Temperature of the inner-core boundary of the earth: Melting of iron at high pressure from first-principles coexistence simulations. Phys. Rev. B, 79:060101, 2009.

D. Alf. Phon: A program to calculate phonons using the small displacement method. Comp. Phys.

Comm., 180:2622Ц2633, 2009.

T. W. Barbee, L. Seaman, R. Crewdson, and D. R. Curran. Dynamic fracture criteria for ductile and brittle metals. J. of Materials, 7:393Ц401, 1972.

A. B. Belonoshko, R. Ahuja, and B. Johansson. Quasi ab initio molecular dynamic study of ferrum melting. Phys. Rev. Lett, 84:3638, 2000.

P. Brommer and F. Gahler. Potfit: effective potentials from ab initio data. Modelling Simulation Mater. Sci. Eng., 15:295, 2007.

J. Chen, D. Tzou, and J. Beraun. Appl. Phys. Lett., 46:307Ц316, 2006.

B. J. Demaske, V. V. Zhakhovsky, N. A. Inogamov, and I. I. Oleynik. Ablation and spallation of gold films irradiated by ultrashort laser pulses. Phys. Rev. B, 82(6):064113, 2010.

J. W. Edwards, R. Speiser, and H. L. Johnston. High temperature structure and thermal expansion of some metals as determined by x-ray diffraction data. i. platinum, tantalum, niobium, and molybdenum. J. App. Physics, 22:424, 1951.

R. Ernstorfer, M. Hard, C. Hebeisen, G. Sciaini, T. Dartigalongue, and R. D. Miller. The formation of warm dense matter: Experimental evidence for electronic bond hardening in gold. Science, 323:

1033, 2009.

R. A. Evarestov and M. V. Losev. All-electron lcao calculations of the lif crystal phonon spectrum:

Influence of the basis set, the exchange-correlation functional, and the supercell size. J. Comp.

Chem., 30:2645Ц2655, 2009.

A. Y. Faenov, N. A. Inogamov, V. V. Zhakhovskii, V. A. Khokhlov, K. Nishihara, Y. Kato, M. Tanaka, T. A. Pikuz, M. Kishimoto, M. Ishino, M. Nishikino, T. Nakamura, Y. Fukuda, S. V. Bulanov, and T. Kawachi. Low-threshold ablation of dielectrics irradiated by picosecond soft x-ray laser pulses.

Applied Physics Letters, 94(23):231107, 2009.

Y. Gan and J. Chen. Integrated continuum-atomistic modeling of nonthermal ablation of gold nanofilms by femtosecond lasers. Appl. Phys. Lett., 94:201116, 2009.

G. V. Garkushin, S. V. Razorenov, and G. I. Kanel. Submicrosecond strength of the d16t aluminum alloy at room and elevated temperatures. Physics of the solid state, 50:805, 2008.

J. A. Gorman, D. S. Wood, and J. Vreeland, T. Mobility of dislocations in aluminum. J. Appl. Phys., 40(2):833Ц841, 1969.

S. Han, L. A. Zepeda-Ruiz, G. J. Ackland, R. Car, and D. J. Srolovitz. Self-interstitials in vanadium and molybdenum. Phys. Rev. B, 66:220101(R), 2002.

A. Hikata, R. A. Johnson, and C. Elbaum. Interaction of dislocations with electrons and with phonons.

Phys. Rev. B, 2:4856, 1970.

W. Hu, Y. C. Shin, and G. King. Energy transpor tanalysis in ultrashort pulse laser ablation through combined molecular dynamics and monte carlo simulation. Phys. Rev. B, 82:094111, 2010.

N. A. Inogamov, V. V. Zhakhovskii, S. I. Ashitkov, Y. V. Petrov, M. B. Agranat, S. I. Anisimov, K. Nishikhara, and V. E. Fortov. Nanospallation induced by an ultrashort laser pulse. Journal of Experimental and Theoretical Physics, 107:1, 2008.

D. S. Ivanov and L. V. Zhigilei. Combined atomistic-continuum modeling of short-pulse laser melting and disintegration of metal films. Phys. Rev. B, 68:064114, 2003.

K. A. Jackson. Kinetic Processes: Crystal Growth, Diffusion, and Phase Transitions in Materials.

Wiley-VCH, Weinheim, 2004.

G. I. Kanel, S. V. Razorenov, K. Baumung, and J. Singer. Dynamic yield and tensile strength of aluminum single crystals at temperatures up to the melting point. J. Appl. Phys., 90:136Ц143, 2001.

G. I. Kanel, S. V. Razorenov, and V. E. Fortov. Shock-wave compression and tension of solids at elevated temperatures: superheated crystal states, pre-melting, and anomalous growth of the yield strength. J. Phys.: Condenced Matter, 16:S1007ЦS1016, 2004.

G. Kresse and J. Furthmuller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 54:11169, 1996.

A. Laio, S. Bernard, G. L. Chiarotti, S. Scandolo, and E. Tosatti. Physics of iron at earthТs core conditions. Science, 287:1027, 2000.

Z. Lin, L. V. Zhigilei, and V. Celli. Electron-phonon coupling and electron heat capacity of metals under conditions of strong electron-phonon nonequilibrium. Phys. Rev. B, 77:075133, 2008.

X.-Y. Liu, F. Ercolessi, and J. B. Adams. Aluminium interatomic potential from density functional theory calculations with improved stacking fault energy. Modelling Simul. Mater. Sci. Eng., 12:665, 2004.

J. F. Lutsko, D. Wolf, S. R. Phillpot, and S. Yip. Molecular-dynamics study of lattice-defect-nucleated melting in metals using an embedded-atom-method potential. Phys. Rev. B, 40(5):2841Ц2855, Aug 1989. doi: 10.1103/PhysRevB.40.2841.

N. D. Mermin. Thermal properties of inhomogeneous electron gas. Phys. Rev., 137:A1441ЦA1443, 1965.

Y. Mishin, D. Farkas, M. J. Mehl, and D. A. Papaconstantopoulos. Interatomic potentials for monoatomic metals from experimental data and ab initio calculations. Phys. Rev. B, 59(5):3393 - 3407, 1999.

Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, A. F. Voter, and J. D. Kress. Structural stability and lattice defects in copper: Ab initio, tight-binding, and embedded-atom calculations. Phys. Rev.

B., 63:224106, 2001.

D. Nguyen-Manh, A. P. Horsfield, and S. L. Dudarev. Self-interstitial atom defects in bcc transition metals: Group-specific trends. Phys. Rev. B, 73:020101, 2006.

S. J. Plimpton. Fast parallel algorithms for short-range molecular dynamics. J. Comp. Phys., 117:

1Ц19, 1995.

S. G. Psakhie, K. P. Zolnikov, and D. S. Kryzhevich. Elementary atomistic mechanism of crystal plasticity. Physics Letters A, 367(3):250Ц253, 2007.

V. Recoules, J. Clrouin, G. Zrah, P. M. Anglade, and S. Mazevet. Effect of intense laser irradiation on the lattice stability of semiconductors and metals. Phys. Rev. Lett., 96:55503, 2006.

C. Schafer and H. M. Urbassek. Metal ablation by picosecond laser pulses: A hybrid simulation. Phys.

Rev. B, 66:115404, 2002.

A. Y. Vorobyev and C. Guo. Enhanced absorptance of gold following multipulse femtosecond laser ablation. Phys. Rev. B, 72:195422, 2005.

Авторефераты по всем темам  >>  Авторефераты по физике