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

2
НЕЛИНЕЙНЫЙ ОТКЛИК МАТЕРИАЛОВ НА МИКРОМАСШТАБНОМ УРОВНЕ ПРИ ВЫСОКОЭ11ЕРГЕТИЧЕСКИХ ВОЗДЕЙСТВИЯХ
Введение 4
1. Методология численного исследования закономерностей отклика материала на атомном уровне при динамическом нагружении 13
1.1 Особенности применения метода молекулярной динамики для изучения нелинейных эффектов
1.2 Парное приближение при описании межатомного взаимодействия
1.3 Многочастичные потенциалы межатомного взаимодействия
1.4 Специфика задания начальных и граничных условий в методе молекулярной динамики
2. Нелинейный отклик и генерация динамических дефектов в области границ раздела при динамическом нагружении и в пост-нагруженном
материале 27
2.1 Особенности нелинейного отклика материалов и генерации дефектов в условиях одноосного растяжения 27
2.2 Генерация вихревого коллективного движения атомов
в зернограничной области при сдвиговом воздействии 44
2.3 О возможности высокоскоростного перемещения границ зерен
при динамическом нагружении 51
2.4 Релаксационные процессы в пост нагруженном материале 60
2.5 Построение диаграмм состояния типа температура-атомная плотность в металлах 68
2.6 Влияние изменений атомной плотности на картину фазовых равновесий в бинарных сплавах 88
13
17
20
24
3
3. Нелинейный отклик материала на микроуровне при высокоэнергетических импульсных воздействиях
106
3.1 Уединенные волны в одномерной цепочке атомов с нелинейными взаимодействиями
3.2 Генерация нелинейных уединенных импульсов в 30 кристаллах при высокоскоростном механическом нагружении
3.3 Формирование уединенных импульсов при ионной имплантации. Эффект дальнодействия
3.4 Нелинейные эффекты при высокоскоростном термическом воздействии
4. Взаимодействие нелинейных уединенных импульсов с дефектами структуры
4.1 Генерация “тепловых пятен” в материале при взаимодействии уединенных импульсов с комплексами дефектов
4.2 Особенности прохождения уединенных импульсов через области с протяженными дефектами структуры
4.3 Эффект отрыва нано-фрагментов материала от свободной поверхности уединенными импульсами
4.4 Влияние протяженных дефектов структуры на характер взаимодействия уединенных импульсов со свободной поверхностью
4.5 Нелинейный перенос энергии при локальном воздействии на материал. Перспективная технология формирования наноскопических многослойных покрытий
Основные результаты и выводы
106
148
163
166
178
178
195
210
223
234
246
Литература
249
ВВЕДЕНИЕ
Компьютерное исследование процессов, происходящих в материалах в условиях экстремальных внешних воздействий, имеет важное научно-практическое значение, поскольку ведет к углублению знаний о закономерностях поведения материала в нелинейной области и позволяет более эффективно решать сложные задачи, связанные с компьютерным конструированием материалов [1-3]. Первоначально интерес к изучению поведения материалов в условиях высокоэнергетических воздействий ориентировался на прогнозирование прочности материалов при высокоскоростных ударах, облучениях различными корпускулярными потоками и т.д. В последние годы в связи с развитием представлений физической мезомеханики о поведении материалов как сложной иерархически организованной системы все более актуальным становится изучение высокоскоростной деформации и разрушения на микроскопическом (атомном) уровне, что важно для более глубокого понимания данных процессов на разных масштабных уровнях [4-6]. В силу названных причин компьютерное моделирование поведения материалов под воздействием таких внешних факторов, как высокоскоростное температурное и механическое нагружения, радиационное облучение, воздействие высокоэнергетическими электронными пучками и т.д., получает в последние годы всё более широкое распространение [1-23].
Как правило, при решении подобных задач на микроскопическом (атомном) уровне в материаловедении используется метод молекулярной динамики. Данный метод основан на решении уравнений движения частиц (атомов), взаимодействие между которыми задается по определенному закону. Спектр приложений этого метода к решению научных задач достаточно широк. Метод молекулярной динамики может быть успешно применен к изучению поведения материалов в условиях высокоэнергетических воздействий.
Средства компьютерного моделирования, основанные на приложении метода молекулярной динамики, позволяют детально исследовать атомные
5
механизмы, отвечающие за проявление тех или иных закономерностей материала в условиях нагружения, и проследить за изменением его свойств и отклика в динамике. Учитывая микроскопический размер исследуемых объектов (десятки ангстрем), короткие временные интервалы наблюдения (наносекунды или даже доли наносекунд), а также сложную систему внешнего нагружения, результаты компьютерного моделирования являются практически единственным источником детальной информации о поведении материала на атомном уровне. Прежде всего, это относится к изучению нелинейных эффектов, которые, как правило, сопровождают или даже определяют процессы, инициированные высокоэнергетическими воздействиями. Методы моделирования позволяют выявлять и анализировать вихревой характер поведения атомной системы в сложных условиях нагружения гетерогенных и композитных материалов, особенно в зонах несоразмерной деформации [1-3].
Таким образом, актуальность исследований, проведенных в настоящей работе, связана с получением новых данных о поведении материалов при высокоэнергетическом нагружении на микромасштабном уровне. Эти данные важны как для более глубокого понимания нелинейных процессов в конденсированных средах, так и для их использования в прикладных аспектах науки о материалах.
Цель работы состояла в изучении нелинейного отклика ЗЭ материалов в условиях высокоэнергетических воздействий на атомном уровне. Для достижения данной цели были сформулированы следующие задачи:
1) изучить характерные особенности вихревого отклика атомной системы в гетерогенных материалах при высокоскоростном нагружении;
2) провести моделирование поведения материала в зоне несовместной деформации при динамическом нагружении;
3) исследовать возможность генерации локализованных нелинейных возмущений в трехмерных кристаллических материалах при импульсном высокоэнергетическом механическом и термическом воздействиях;
4) исследовать взаимодействие нелинейных возмущений с дефектами структуры материала;
6
5) провести моделирование взаимодействия нелинейных возмущений, сформированных высокоэнергетическим импульсным нагружением, со свободной поверхностью материала.
Научная новизна. В работе получены следующие новые результаты:
- показано, что вихревая мода в поведении атомной системы может играть важную роль при динамическом нагружении кристаллитов, обеспечивая, в частности, такие важные процессы аккомодации, как аномально высокая скорость перемещения границ зерен и фрагментация;
впервые показано, что в трёхмерном кристаллите при высокоэнергетическом импульсном воздействии могут формироваться нелинейные уединенные импульсы. Изучен характер их взаимодействия друг с другом и дефектами структуры;
- обнаружена возможность генерации "тепловых пятен” при прохождении нелинейных уединенных импульсов через вакансионные комплексы и исследованы их параметры;
- впервые показана возможность процессов отрыва при выходе нелинейных уединенных импульсов на свободную поверхность. Изучены параметры процесса отрыва (скорости отлета фрагментов, их размеры и т.д.) в зависимости от количества уединенных импульсов, их амплитуд, кристаллографической ориентации свободной поверхности, наличия протяженных дефектов типа границ зерен в приповерхностной области и т.д.;
- предложены и обоснованы новые подходы по практическому использованию нелинейных уединенных импульсов.
Научная и практическая ценность.
Проведенные исследования расширяют представления о возможностях нелинейного отклика материалов на микроуровне в условиях высокоэнергетических воздействий.
На основе полученных результатов описаны атомные механизмы фрагментации материалов, перемещения межзеренной границы, локального температурного разогрева, структурных перестроек в области границ зерен, происходящих под воздействием динамического нагружения.
7
Результаты моделирования взаимодействия нелинейных возмущений с комплексами дефектов могут быть использованы при анализе возможных механизмов твердофазных химических реакций и процессов механического перемешивания компонентов материала при высокоскоростном воздействии.
Закономерности распространения уединенных импульсов в реальных материалах представляют интерес при разработке новых методов неразрушающего контроля накопления дефектов в эксплуатируемых материалах.
Эффекты отрыва при взаимодействии нелинейных уединенных импульсов со свободной поверхностью могут быть положены в основу нового подхода по нанесению сверхтонких многослойных покрытий.
Положения, выносимые на защиту:
1. Вихревой характер поведения атомной системы кристаллита, посредством которого может осуществляться аномально высокая скорость перемещения границ зерен в условиях высокоэнергетических воздействий.
2. Особенности структурных изменений и атомные механизмы аккомодации на микроуровне в зонах несовместной деформации при динамическом нагружении.
3. Генерация и свойства нелинейных уединенных импульсов в 3D кристаллах при импульсном высокоэнергетическом нагружении.
4. Особенности взаимодействия нелинейных уединенных импульсов с дефектами кристаллической структуры.
5. Процессы отрыва фрагментов материала при выходе нелинейных уединенных импульсов на свободную поверхность и метод получения сверхтонких многослойных покрытий.
Апробация работы. Материалы диссертации докладывались на: XI International Conference HERF (Ljublijana, Yugoslavia, 1989); 4-ой международной конференции "Компьютерное конструирование перспективных материалов и технологий” CADAMT-95, (Томск, Россия, 1995); Second Sino-Russian Symposium “Advanced Materials and Processes”
8
(Xian, China, 1995); 8-ой международной конференции “Уравнения состояния” (Эльбрус, Россия, 1992); международной конференции "High Power Lasers - Science and Engineering" (Карловы Вары, Чехия, 1995); международной конференции "Metallurgical and Materials Applications of Shock-Wave and High-Strain-Rate Phenomena" (El Paso, Texas, USA, EXPLOMET 1995); международной конференции “Computer-Aided Design of Advanced Materials and Technologies — CADAMT” (Байкальск, 1997); межгосударственной конференции "Новые численные методы упругости и пластичности" (Новосибирск, 1997); International Conference “Movable cellular automata method: Foundation and Application” (Ljublijana, Slovenia, 1997); 11th Conference on Radiation Physics and Chemistry of Condensed Matter (Tomsk, Russia, 2000); 5th Conference on Modification of Materials with Particle Beams and Plasma Flows (Tomsk, Russia, 2000); 6th International Scientific and Practical Conference of Students, Post-Graduates and Young Scientists “Modern Techniques and Technology” (Tomsk, Russia, 2000); Euromech Colloquium “Fracture Aspects in Manufacturing" (Moscow, Russia 2000); VI International Conference “Computer-Aided Design of Materials and Technologies" (Tomsk, Russia, 2001).
Основные результаты диссертации опубликованы в 36 работах, перечень их наименований представлен в списке цитируемой литературы.
Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения и списка цитируемой литературы, содержит 104 рисунка, 2 таблицы, библиографический список из 209 наименований -всего 267 страниц.
Во введении обоснована актуальность исследуемой проблемы, сформулирована цель работы, перечислены новые результаты, раскрыта их практическая ценность, представлены положения, выносимые на защиту, описана структура диссертации.
В первой главе диссертации достаточно подробно излагается методология компьютерного моделирования отклика материала на
9
микроскопическом уровне в условиях динамического нагружения. Для решения подобного рода задач наиболее подходящим является метод молекулярной динамики, который базируется на хорошо определенном микроскопическом описании физической системы. Данная система, как правило, состоит из многих тел и описывается с помощью лагранжиана, гамильтониана или уравнениями движения Ньютона. Расчет свойств физической системы, особенностей ее поведения в условиях разнообразных внешних воздействий проводится на основе компьютерного моделирования. При этом ограничения компьютерных экспериментов связаны с ограничениями на количество рассматриваемых частиц в моделируемой системе и ограничениями на продолжительность времени, в течение которого изучается поведение системы частиц. Данные ограничения вызваны как быстродействием, так и памятью существующих компьютеров. В настоящей работе в качестве изучаемой системы будет выбрано твердое тело, имеющее кристаллическую структуру. Важную роль при молекулярно-динамических расчетах играет тип используемых потенциалов межатомного взаимодействия. В главе приведены выражения для эмпирических, парных и многочастичных потенциалов межатомного взаимодействия и определены области их применимости, отмечены особенности их применения при решении тех или иных задач. Для корректного решения системы дифференциальных уравнений в методе молекулярной динамики необходимо правильно смоделировать начальные условия и выбрать соответствующие граничные условия. Эти вопросы обсуждаются в последнем параграфе первой главы.
Во второй главе исследуется поведение гетерогенных материалов при одноосном высокоскоростном растяжении. Хотя этой теме посвящено большое число экспериментальных и теоретических работ, тем не менее особенности отклика материала на атомном уровне еще недостаточно выяснены. В данной главе описаны результаты моделирования нелинейного отклика 30 материалов при динамическом нагружении вдоль различных кристаллографических направлений со скоростями, лежащими в интервале от 1 до 100м/с, изучено влияние концентраторов напряжения
10
(затравочных микротрещин) на особенности поведения материала в условиях одноосного динамического растяжения. Показано, что
аккомодационные процессы в материале в условиях высокоскоростного нагружения являются инерционными и продолжают активно протекать после снятия внешнего воздействия.
Проанализированы результаты моделирования нелинейного
поведения материалов, содержащих межзеренные границы, на высокоскоростное сдвиговое нагружение. Показано, что одним из основных способов отклика материала на такое нагружение является формирование вихревого движения атомной системы в зернограничной области. Описаны размеры вихрей, условия их генерации, показано, что вихревое движение атомных групп может приводить к высокоскоростному перемещению границ зерен в материале. Исследованы структурные перестройки в
зернограничной области как в процессе нагружения, так и после прекращения внешнего воздействия. Предложена методика построения новых фазовых диаграмм состояния для равновесных состояний, у которых одной из переменных является атомная плотность. Данная методика проиллюстрирована расчетами диаграмм состояния температура-атомная плотность для металлов и трехмерными диаграммами типа температура-атомная плотность-концентрация для бинарных металлических систем. Изучено влияние всестороннего давления на изменение диаграмм состояния температура-атомная плотность для бинарных систем, характеризующихся неограниченной взаимной растворимостью в твердой фазе.
Третья глава посвящена нелинейным уединенным импульсам, которые генерируются в материалах в результате высокоэнергетического воздействия. Изучены особенности генерации и свойства уединенных импульсов как в одномерных цепочках, так и в ЗИ материалах. Показано, что по своим свойствам уединенные импульсы схожи с солитонами. Скорость распространения уединенных импульсов зависит от их амплитуды, после взаимодействия друг с другом они восстанавливают свою форму, в материалах с идеальной кристаллической структурой они
11
переносят заключенную в них энергию в баллистическом режиме. Уединенные импульсы могут генерироваться как в результате высокоэнергетического механического нагружения, так и в результате высокоскоростного разогрева или ионной бомбардировки свободной поверхности. Способность уединенных импульсов в баллистическом режиме переносить энергию позволяет глубже понять природу “эффекта дальнодействия”, наблюдаемого при ионной имплантации металлических материалов.
В четвертой главе исследуются особенности взаимодействия нелинейных уединенных импульсов с дефектами структуры: комплексами вакансий, областями с пониженной атомной плотностью, границами зерен. Показано, что в результате взаимодействия с вакансионным комплексом в материале возможна генерация “теплового пятна”. Изучены характеристики генерируемых тепловых пятен (их размеры, изменение температуры со временем), а также показана возможность их генерации как результат взаимодействия с близкорасположенным неравновесным дефектом. Исследованы особенности прохождения уединенных импульсов через границы зерен, проведена оценка относительного изменения энергии уединенных импульсов после взаимодействия с зернограничной областью.
Показано, что выход на свободную поверхность уединенных импульсов с амплитудами выше порогового значения приводит к отрыву приповерхностных фрагментов. Геометрия отрыва зависит как от кристаллографического направления распространения импульсов, так и от наличия в приповерхностной области дефектов структуры. Скорости отлета оторвавшихся фрагментов могут достигать нескольких километров в секунду. При наличии в приповерхностной области границ раздела существует тенденция к отрыву вдоль границы раздела. Способность к отрыву вдоль границы усиливается, если в материале сформирован пакет уединенных импульсов. Проведено исследование особенностей эффекта отрыва приповерхностных фрагментов при локальном высокоэнергетическом нагружении. Указано, что способность частичного рассеивания уединенных импульсов на дефектах структуры представляет
12
интерес для решения задач по неразрушающему контролю за накоплением повреждений в материале при нагружении. Предложен новый подход по нанесению сверхтонких многослойных покрытий, основанный на эффекте отрыва при взаимодействии уединенных импульсов со свободной поверхностью.
В заключении диссертации приводятся основные результаты и выводы.
Глубокую признательность и благодарность автор выражает научному консультанту профессору Псахье Сергею Григорьевичу за ценный опыт и знания, полученные за годы многолетнего сотрудничества, за постоянное внимание и ценные советы при работе над диссертацией. Автор благодарит академика Панина Виктора Евгеньевича за постоянное внимание и ценное обсуждение результатов диссертационной работы. Автор признателен всем сотрудникам лаборатории компьютерного конструирования материалов ИФПМ СО РАН за плодотворное сотрудничество, ценные обсуждения и помощь в работе. Автор считает своим долгом поблагодарить жену, Зольникову Л.М., внимание и поддержка которой помогали ему во все времена.
ГЛАВА 1. МЕТОДОЛОГИЯ ЧИСЛЕННОГО ИССЛЕДОВАНИЯ ЗАКОНОМЕРНОСТЕЙ ОТКЛИКА МАТЕРИАЛА НА АТОМНОМ УРОВНЕ ПРИ ДИНАМИЧЕСКОМ НАГРУЖЕНИИ
1.1 Особенности применения метода молекулярной динамики для изучения нелинейных эффектов
Аналитическое изучение свойств материалов на основе использования знаний о межатомных потенциалах взаимодействия представляет собой неразрешимую задачу молекулярной теории. Основным препятствием в проведении аналитических расчетов является проблема многих тел. Решением данной проблемы может служить компьютерное моделирование систем, состоящих из большого числа частиц. При этом ограничения компьютерных экспериментов связаны с ограничениями на количество рассматриваемых частиц в моделируемой системе и ограничениями на продолжительность времени, в течение которого изучается поведение системы частиц. Данные ограничения вызваны как быстродействием, так и памятью существующих компьютеров. Одним из основных численных методов в статистической механике, позволяющего рассматривать системы с большим количеством частиц, является метод молекулярной динамики.
В рамках метода молекулярной динамики моделируемый образец рассматривается как система N классических частиц (атомов), взаимодействующих по определенному закону, который задается потенциалом взаимодействия. Состояние такой системы описывается 6Л/-мерным вектором и, образованным пространственными координатами Х(Х = х,у,г) и компонентами скоростей У(У = УяУуу.) всех частиц системы.
Тогда, согласно законам движения Ньютона, эволюция совокупности N атомов во времени описывается системой 6Л/ обыкновенных дифференциальных уравнений первого порядка типа [24]:
14
для функции и = £/(/) с некоторыми начальными условиями г/(/°) = и0.
Интегрируя уравнение (1.1) по малому временному шагу д/, можно
получить связь между состояниями и"*1 и ип в двух соседних точках на временной оси:
Для аппроксимации интеграла по времени в правой части уравнения
(1.2) в большинстве случаев достаточно использования простейшей схемы Эйлера первого порядка [25]. Согласно схеме Эйлера:
Запишем дифференциальные уравнения, определяющие эволюцию во времени, выбранной совокупности классических частиц. Для ьй частицы они имеют вид:
где V,х.,т, - скорость, пространственная координата и масса ьй частицы;
(рь - потенциал взаимодействия ьй частицы с]-ой; г1{ = = |/~, .
Чтобы избежать разностного дифференцирования по пространству, при вычислениях целесообразно использовать не потенциалы, а скалярную силу:
дг
тогда система (1.4) запишется в виде:
(1-2)
и'*1 =и" -/(£/",/")-Д/
(1.3)
л
<1У, _ 1 ^ <5, ’
си т, дг,г Ж-д
(1.4)
Л
(1.5)
Согласно (1.5) разностная схема запишется в виде:
Таким образом, появляется возможность рассчитывать состояние системы для каждого из дискретного множества значений времени п-1,2,3... Это позволяет проводить исследование поэтапной эволюции выбранной системы во времени. Для обеспечения сохранения энергии полезно использовать третий закон Ньютона таким образом, чтобы изменение импульса /-й частицы, обусловленное взаимодействием с у-ой, было равно изменению импульса с обратным знаком )-й частицы, обусловленному взаимодействием с /-ой частицей.
Необходимо отметить, что при микроскопических расчетах обычно используется атомная система единиц [26,27], в которой масса, заряд электрона, боровский радиус и постоянная Планка равны единице.
При решении разностных уравнений (1.6) основным параметром является временной шаг . Приращение Аг следует выбирать с условием, чтобы для всех '! (1<='к=п) выполнялось соотношение:
Л/ <КЩ'
где а - характерное расстояние взаимодействия.
Метод молекулярной динамики нашел применение для решения колоссально широкого круга задач. Среди них можно упомянуть моделирование фазовых переходов твердое тело - жидкость и жидкость -твердое тело [28-30], исследования вязкого течения жидкости [31-33], поведения материала на микроскопическом уровне при ударно-волновом и радиационном нагружении материалов [34-44] и т.д..
Из (1.6) следует, что поведение атомной подсистемы может быть описано, если известен потенциал межатомного взаимодействия. Простейший потенциал взаимодействия определяется в модели твердых сфер [45,46]. Более реалистичным подходом является использование непрерывной потенциальной функции, в которой ее первая и вторая производные непрерывны, как, например, в [47-52]. К таким потенциалам можно отнести потенциал Борна-Майера:
16
<р(г) = Аехр(-В/г),
(1.7)
где А и В - подгоночные параметры, г - расстояние между взаимодействующими частицами. Потенциалы (1.7) являются "отталкивательными" и пригодны для описания системы "мягких" сфер.
Потенциалы, описывающие отталкивание и притяжение между взаимодействующими частицами, были предложены Леннардом-Джонсом [41]:
где параметр е - глубина потенциальной ямы, сг - определяет значение равновесного расстояния; и Морзе [52]:
где А и а параметры моделируемого материала.
Отметим, что потенциалы Леннарда-Джонса наиболее часто применяются для описания поведения инертных газов в различных агрегатных состояниях, а потенциалы Морзе используются для моделирования поведения металлов и полупроводников. В последующем было предложено множество различных эмпирических потенциалов, представляющих собой полиномы, а также имеющие экспоненциальную зависимость от межатомного расстояния. Параметры полуэмпирических и эмпирических потенциалов часто подгоняют по свойствам, зависящим только от объема.
Существенный прогресс в моделировании поведения материалов на атомном уровне начался с развитием псевдопотенциальной теории и расчетов эффективных потенциалов на ее основе. Значение этой теории для компьютерных экспериментов трудно переоценить, поскольку корректность результатов моделирования, проведенных в рамках метода молекулярной динамики, во многом определяется потенциалами межчастичного взаимодействия. Поэтому в следующих параграфах рассмотрим подробнее те потенциалы, которые будут использованы в численных расчетах данной диссертации.
(1.8)
<р(г) = л|ехр|- 2а|г - ))- 2ехр(- а[г - |,
(19)
17
1.2. Парное приближение при описании межатомного взаимодействия
С конца шестидесятых годов прошлого столетия в физике конденсированного состояния широкое применение получил метод псевдопотенциала. Во многом это связано с простотой самого метода, позволяющего описать с достаточно высокой степенью точности многие свойства металлов и сплавов, находящихся в различных агрегатных состояниях.
Теория псевдопотенциала основана на следующих приближениях [26,53]:
1. Адиабатическое приближение. Суть его состоит в пренебрежении движением ионной системы. Поскольку кинетическая энергия ионов и электронов проводимости примерно равны, то ввиду большой массы ионов их скорости движения много меньше скорости движения электронов.
2. Приближение самосогласованного поля. В этом приближении прямое взаимодействие электронов заменяется некоторым усредненным полем. Это усредненное поле оказывает влияние на движение каждого электрона и, в свою очередь, само зависит от движения каждого электрона. Поэтому данное поле нужно вычислять самосогласованным образом.
3. В основе третьего приближения лежит подразделение всех электронов на электроны "остова" и электроны зоны проводимости. Предполагается, что состояния электронов "остовов" являются локализованными, совпадающими с соответствующими электронными состояниями изолированных атомов. В то же время состояния электронов зоны проводимости считаются нелокализованными.
Важной чертой метода псевдопотенциала является возможность введения параметра малости (отношение экранированного псевдопотенциала "остова" мало по сравнению с энергией Ферми). Это, в свою очередь позволяет применить в полном объеме теорию возмущений к электронам проводимости. Во втором порядке малости теории возмущений полная энергия моделируемой системы может быть представлена [26,53]:
18
'Г'і
О
(1.10)
где /'’(Пф) - вклад, зависящий только от объема о0; второе слагаемое <р„
определяет вклад в полную энергию системы, зависящий от взаимного расположения атомов.
Эффективный потенциал парного межатомного взаимодействия определяется следующим образом:
% = —1+— I Кля)- ж я Ж/.
г 2/Т' аг
(1.11)
0 ЯГ
где - валентности I - ого и ] - ого ионов, соответственно, ц - вектор в обратном пространстве, ^{ц) - характеристическая функция зонной
структуры, зависящая только от природы химического элемента и не зависящая от кристаллической структуры:
/’ЬМ-ЬчМаМ«). (1.12)
где - псевдопотенциалы ионов, £(</)- диэлектрическая функция,
которая учитывает обменно-корреляционные эффекты в однородном электронном газе. Выражение для с (</) может быть записано в виде [26]:
* (?)881 “ 7ГТІ1 - ЛяЪс (я)
В выражении (1.12) и (1.13) х(ч) имеет следующий вид:
9 +
(113)
*(?)=-
/ \ /
32
[2^7
\ > /
І+і£Ьі1щ
2 8 дк,
Я ~ 2к}
(1.14)
где к/ =(зж2г/П0)п - радиус сферы Ферми, Е, - энергия Ферми, Е, - к) 12.
Функция /(({) в выражении (1.13) учитывает поправки к кулоновскому взаимодействию электронов за счет эффектов обмена и корреляции в электроном газе и ведет к уменьшению этого взаимодействия. Точного выражения для функции /(</) в многоэлектронной теории не удается получить, поэтому используется ряд приближенных формул для этой функции [54].
19
В настоящей работе для учета эффектов обмена и корреляции в диэлектрической функции использовалось приближение Гелдарта-Воско [55]:
В качестве псевдопотенциала, описывающего рассеяние электронов проводимости на ионах, выбирается функция с некоторым числом параметров, которая строится по определенным правилам. Эта функция называется модельным псевдопотенциалом. В ряде расчетов, выполненных в настоящей работе, в качестве модельного локального псевдопотенциала использовался двухпараметрический потенциал Краско-Гурского [56], Фурье-представление которого имеет вид:
где а,гс - подгоночные параметры.
Параметры псевдопотенциала, как и в работах [57], определялись из условия наилучшего совпадения со следующими экспериментальными данными:
экспериментальное значение; р - давление; п - равновесный объем, приходящийся на один атом, при нулевых значениях температуры и давления.
(1.15)
(1.16)
р(О) = 0
рассчитанное значение модуля упругости;
20
Выбор модельного потенциала в форме Краско-Гурского и расчета параметров потенциала согласно условию (1.16) позволяет хорошо описывать статические и динамические свойства простых металлов.
1.3. Многочастичные потенциалы межатомного взаимодействия
Несмотря на то, что потенциалы межатомного взаимодействия, полученные в рамках псевдопотенциальной теории, позволяют описывать с достаточно хорошей точностью широкий круг свойств металлов и сплавов, в ряде случаев их применение становится некорректным; в частности, при моделировании протяженных дефектов структуры, внутренних границ раздела типа межфазных и межзеренных границ, свободной поверхности и т.д.. Следует отметить, что атомная и электронная структура свободной поверхности в принципе не описывается в рамках модельного парного псевдопотенциала, поскольку, как показали расчеты [54], первая координационная сфера для них располагается на отталкивательной ветви. Поэтому эффект поджатия поверхностных атомных плоскостей, который характерен для большинства кристаллографических направлений, не может быть описан в рамках парного приближения для потенциалов межатомного взаимодействия. В то же время экспериментальные и теоретические исследования последних десятилетий убедительно показали, что физико-химические процессы, происходящие на поверхности и в области границ раздела, играют важнейшую роль в зарождении дефектов структуры, зарождении и развитии разрушения, роста и формообразовании кристаллов и т.д.. Кроме того, поверхности материалов во многом определяют процессы спекания порошков, коррозию, окисление, трение и т.д.. Наиболее подходящими для моделирования поверхностных явлений в настоящее время являются самосогласованные расчеты из первых принципов. К таким расчетам можно отнести построение многочастичных потенциалов взаимодействия, выполненное в рамках метода погруженного атома [58,59]. Моделирование поверхностных явлений с использованием многочастичных потенциалов имеет важное
21
значение еще и потому, что в экспериментальных измерениях поверхностной энергии наблюдается значительный разброс [2,60].
В методе погруженного атома [58,59] каждый атом представляется помещенным в почти однородный электронный газ, а энергетическая добавка, связанная с изменением энергии при добавлении атома в систему, входит как составная часть в выражение для полной энергии системы. В [58] Дау показал, что плотность неоднородного электронного газа в этом случае должна не сильно отличаться от суперпозиции атомных плотностей. Это обусловило перспективность применения метода погруженного атома для расчета потенциалов межатомного взаимодействия в металлах и их сплавах. Данный метод достаточно хорошо работает для металлических систем с электронной плотностью, незначительно отклоняющейся от суперпозиции атомных плотностей, в частности, для переходных металлов с почти заполненной с! - оболочкой (N1, Си, Рс1, Ад, Р\, Аи).
В общем случае полная энергия кристалла в методе погруженного атома может быть записана следующим образом [58,59]:
Е = \ТфММ1г(р.) (117)
где Ф;(д,)- парные потенциалы межатомного взаимодействия: -
расстояние между положениями 1-ого и ]-ого атома; Р(р,)~ функция погружения, т.е. изменение энергии кристалла при помещении в него /-го атома, ад- плотность, приходящая на ядро /-го атома со стороны остальных атомов системы:
Д = 2>ф(-Я, |), (1.18)
./
где р*(г)- парциальные атомные плотности. Представление электронной плотности в виде суперпозиции атомных плотностей является приближением. При переходе от изолированного атома к металлу энергия каждого атома меняется незначительно, поэтому возмущающий потенциал, который действует на электроны рассматриваемого атома со стороны окружающих атомов, мал. Соответственно мала и погрешность, возникающая при аппроксимации электронной плотности в кристалле
22
суперпозицией атомных плотностей. В настоящее время в рамках метода погруженного атома применяют различные выражения для представления парных потенциалов межатомного взаимодействия фД/?„) и парциальной
электронной плотности /?"'(/). Так в работе [58,59] парные потенциалы определялись эмпирическими зависимостями вида:
гтглю
, (1.19)
и Я
где Z(R) - эффективный заряд, монотонно уменьшающийся с увеличением расстояния:
г(Я) = г{)(\ + /зг )ехР (-ой) , (1.20)
где 70 - число внешних электронов; у, а, р - параметры.
Значение V бралось равным 1 для №,Рб,Р1 и 2 для Си,Ад,Аи, что позволяло хорошо описывать упругие постоянные [58].
Изменение энергии системы при погружении атомов в электронный газ представлено в виде аддитивных функций погружения Р(р). Функция
погружения описывает взаимодействие между «погруженным» атомом и остальными атомами кристалла через электронную плотность в области «погружения» атома. Одним из способов нахождения функции погружения является расчет энергии идеального кристалла при его всестороннем сжатии и растяжении, когда в силу симметрии все атомы находятся в
эквивалентных положениях. В [58,61 ] функцию погружения предложено
определять через зависимость величины энергии образования твердого тела от объема, заданную универсальным уравнением состояния Розе-Винетта [62]:
Е, =£0(1 + я*)ехр(-я*), (1.21)
где Е0 - энергия сублимации при равновесном значении параметра решетки а=а0, а* =(а/а^-1)/[£0/(9ЯП)]1/2, а -параметр решетки, О- объем кристалла, приходящийся на один атом, В - модуль всестороннего сжатия.
23
Уравнение состояния (1.21) хорошо описывает зависимость давления от объема при низких температурах для ряда металлов с плотноупакованными кубическими структурами.
В настоящей работе использовались потенциалы, вычисленные в работах [63-65]. В этих работах параметры а и (3 в формуле (1.20), определялись из условия наилучшего описания упругих постоянных и энергий образования вакансий моделируемых металлов.
Отметим, что метод погруженного атома, развитый в [58,59], с успехом применялся для решения широкого круга задач: расчет фононного спектра в объеме твердого тела [66]; исследование жидких металлов [67], дефектов [58,68], структуры поверхности [69-72], поверхностных фононов [73] и термодинамических свойств ГЦК переходных металлов [74] и т.д..
В целях оптимизации компьютерных расчетов в настоящей работе были использованы табличные данные для зависимости потенциала взаимодействия и его производной от квадрата расстояния. Выражение для хк -ой компоненты силы, действующей на атом в ьом узле, в зависимости от квадрата расстояния имеет следующий вид:
24
1.4 Специфика задания начальных и граничных условий в методе молекулярной динамики
Как отмечалось выше, эволюция моделируемой системы определяется уравнениями движения (1.1), а также начальными и граничными условиями. Наиболее распространенными являются варианты жестких и периодических граничных условий. При расчетах с жесткими граничными условиями граничная область представляет собой несколько слоев атомов, жестко фиксированных в узлах кристаллической решетки. Толщина такого слоя определяется эффективным радиусом взаимодействия атомов. Жесткие граничные условия обычно используются для решения статических задач. Следует отметить, что применение жестких граничных условий может приводить к "закачиванию" энергии в моделируемую систему, вследствие взаимодействия подвижных атомов расчетной ячейки с неподвижными атомами граничного слоя. Этот эффект проявляется в большей степени с увеличением размерности моделируемой системы (так как число атомов расчетной ячейки, контактирующих с неподвижным граничным слоем возрастает) и с ростом температуры расчетной ячейки. Влияние граничных эффектов на поведение моделируемой системы увеличивается с ростом расчетного времени. Для минимизации влияния границ на поведение изучаемой системы часто используют периодические граничные условия. Суть периодических граничных условий заключается в том, что любое свойство в моделируемой системе обладает трансляционной симметрией:
А(х) = А(х+1_),
где 1_ - размер расчетной ячейки.
Если в процессе счета координата ной частицы стала больше I (Х|>1_), то ее новая координата определяется как х* = х, - 1_; или, если координата х,<0, то х, = х, + 1_. Для исключения эффектов "самодействия" частиц, эффективный радиус взаимодействия частиц должен быть меньше половины размера моделируемой системы вдоль направления применения периодических граничных условий.