CC BY
16
УДК 621.785.53
О МОДЕЛИРОВАНИИ ДВИЖЕНИЯ КРАЕВЫХ ДИСЛОКАЦИЙ В ОЦК-ЖЕЛЕЗЕ
НАГОРНЫХ И.Л., БУРНЫШЕВ И.Н.
Институт механики Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. В настоящей работе приводится алгоритм оптимизации модели кристалла ОЦК-Fe,
a
содержащего краевую дислокацию — < 111 > {110}, с периодическими граничными условиями и ЕАМ-
потенциалами. Оптимизация позволяет выбрать параметры модели, при которых взаимодействие дислокации с ее образами пренебрежимо мало, а число атомов в моделируемой системе не превышает 50000. Вычисленное значение напряжения Пайерлса-Набарро составляет 75 МПа.
КЛЮЧЕВЫЕ СЛОВА: молекулярная динамика, моделирование, краевая дислокация, вектор Бюргерса, напряжение Пайерлса-Набарро.
ВВЕДЕНИЕ
Подавляющее большинство работ, направленных на исследование дислокационной структуры в металлах, посвящены изучению винтовых дислокаций, а поведение краевых дислокаций в металлах остается малоизученным. Немногие работы, связанные с моделированием краевых дислокаций в металлах, не содержат подробной информации, необходимой для постановки таких вычислительных экспериментов. Такая информация необходима для изучения (моделирования) механизмов разрушения металлов, контактирующих с водородсодержащими средами. Ранее нами методом молекулярной динамики было показано, что водород в зависимости от его концентрации в металле почти на 50 % снижает теоретический предел прочности идеальных кристаллов железа, никеля, алюминия и палладия [1, 2]. В выполненных расчетах не учитывалась дислокационная структура металлов, наличие которой в реальных металлах приводит к снижению предела прочности на несколько порядков.
Целью настоящей работы является формирование молекулярно-динамической модели кристалла, содержащего краевую дислокацию, и обоснование выбора параметров такой модели. Настоящая работа является частью большей работы, посвященной изучению взаимодействия краевых дислокаций с атомами растворенного водорода в металлах. Полные
a
дислокации в ОЦК решетке с вектором Бюргерса — < 111 > (а - параметр решетки) обладают
наименьшей энергией и встречаются чаще остальных. Скольжение в этом случае может осуществляться по плоскостям {110} и {211}. В качестве модельной выбрана дислокация <111>{110}, движущаяся в кристалле ОЦК-железа, не содержащем каких-либо других кристаллических дефектов.
МЕТОДИКА ВЫЧИСЛИТЕЛЬНОГО ЭКСПЕРИМЕНТА
При формировании модели кристалла за основу взята молекулярно-динамическая модель, схематично изображенная на рис. 1. Данная модель представляет собой ориентированный по кристаллографическим направлениям [111], [0-11] и [-211] кристалл, содержащий краевую дислокацию. Периодические граничные условия (111 У) накладываются вдоль линии дислокации ([-211], OZ) и по направлению вектора Бюргерса ([111], OX).
У этой модели есть два недостатка. Первым недостатком такой модели является наличие в системе взаимодействий между дислокацией и её образами по направлению вектора Бюргерса. Такое взаимодействие может привести к ошибочным значениям напряжения Пайерлса-Набарро. Другим недостатком рассматриваемой модели кристалла является наличие в системе свободных поверхностей, искажающих поведение дислокации.
АХ, Дt
в
А
_1_
Н
У[0 1 1]
Ъ\2 1 1] &
ь
X [1 1 1]
Рис. 1. Схема модели кристалла (обозначения приведены в тексте)
Оба этих недостатка могут быть устранены выбором минимально необходимых размеров кристалла вдоль направления вектора Бюргерса ([111], ОХ) и вдоль направления вектора OY [0-11]. Более подробно второй недостаток будет рассмотрен ниже.
На начальном этапе вычислительного эксперимента в кристалле формировалась дислокация. В качестве базовой системы применялся кристалл размерами 10*10x10 нм, ориентированном по кристаллографическим направлениям [111], [0-11], [-211]. Для формирования краевой дислокации 2 < 111 > {110} в идеальном кристалле ОЦК железа удалялись три
полуплоскости, образуя тем самым три экстраплоскости будущей краевой дислокации (рис. 2) [3]. Релаксация такого кристалла под действием внутренних напряжений приводит к
образованию полной краевой дислокации а < 111 >{110}. На втором этапе производили
присвоение нулевых значений скоростям атомов областей В и С, приведенных на рис. 1.
Рис. 2. Проекция координат атомов в области краевой дислокации после релаксации
В дальнейшем значения скоростей атомов области С остаются нулевыми за все время расчета. Затем осуществляли мгновенное смещение атомов области В на величину АХ вдоль оси ОХ с последующей релаксацией системы в течение времени А! Во время релаксации скорости атомов области В остаются нулевыми. Последняя операция повторяется на протяжении всего расчета. Соотношение величин АХ и А1 определяет, прежде всего,
релаксацию кристалла между смещениями атомов области А, при которой не происходит увеличение сдвигового напряжения на плоскостях скольжения сверх напряжения Пайерлса-Набарро.
Все расчеты проводились по схеме Верле с шагом по времени 10 пс для температуры ~ 0,01 K. Потенциалы межатомного взаимодействия рассчитаны в приближении метода погруженного атома (ЕАМ - embedded atom method) [4, 5]. В процессе расчета проводилось измерение сдвигового напряжения 5xy в области, обозначенной пунктиром на рис. 2, по классической теореме вириала и фиксировалось положение дислокации D. Положение дислокации D определялось по координатам атома области А, обладающего наибольшей потенциальной энергией. Расчет останавливался в момент времени, когда дислокация, преодолев расстояние AD = L, оказывалась в положении, соответствующем начальным условиям.
При оптимизации модели в качестве оптимизируемых параметров выбраны: линейный размер кристалла L, определяющий взаимодействие между дислокацией и ее образом по направлению вектора Бюргерса; линейный размер кристалла H, определяющий влияние свободных поверхностей на поведение дислокации; величины AX и At, определяющие релаксацию кристалла после внесения возмущения в систему.
Результаты расчета для системы с параметрами L = 10 нм, H = 10 нм, AX = 0,001 нм, At = 5 пс представлены на рис. 3. На графиках рис. 3 можно определить напряжение Пайерлса-Набарро оп-н. Скачки напряжения сдвига oxy возникают из-за релаксации напряженного состояния, связанного с движением дислокации, что отчетливо видно из рис. 4.
При детальном рассмотрении кривых рис. 3 можно увидеть, что каждый момент возмущения системы отражается на кривой сдвигового напряжения и кривой полной энергии Etot моделируемой системы, что выражается в появлении искажений с шагом 5 пс (рис. 4). Из рис. 3 и 4 видно, что дислокация перемещается дискретно с шагом | b |/3 (b - вектор Бюргерса).
По поведению кривой Etot можно сделать вывод, что релаксация системы после внесения сдвига величиной AX = 0,001 нм происходит за время 1 - 2 пс, релаксация системы после скачка дислокации также происходит в течение 1 - 2 пс. Это означает, что соотношение величины сдвига AX = 0,001 нм и At = 5 пс является достаточным, чтобы система оказалась в устойчивом равновесном состоянии до момента следующего возмущения или перемещения дислокации.
Время, не
Рис. 3. Сдвиговое напряжение и соответствующее ему по времени положение дислокации ДО для модели с параметрами L = 10 нм, Н = 10 нм, АХ = 0,001 нм, At = 5 пс
и
и
И
-52.06955 -
-52.06960
0.23 0.232 0.234 0.236 0.238 0.24 0.242 0.244 0.246 0.248 0.25
Время, нс
Рис. 4. Пример скачка дислокации АБ с сопутствующими изменениями напряжения сдвига ох
и полной энергии системы Еш
Напряжение Пайерлса-Набарро оп-н соответствует минимальному значению напряжения Оху (рис. 4). Из рис. 3 следует, что напряжение Пайерлса-Набарро Оп-н оказывается зависимым от положения дислокации в кристалле и изменяется от минус 140 до минус 40 МПа. Данное явление возникает из-за близкого расположения дислокации к свободной поверхности. Недостаточная удаленность дислокации от поверхности в совокупности с некомпенсированными экстраплоскостями в системе приводит к искривлению плоскостей скольжения {110}.
Обозначенный недостаток удалось устранить увеличением параметра Н в 4 раза. Результаты расчета для системы с параметрами L = 10 нм, Н = 40 нм, АХ = 0,001 нм, Дt = 5 пс представлены на рис. 5. Напряжение Пайерлса-Набарро оп-н составляет 75 МПа. Стоит отметить, что последующее увеличение каждого из линейных размеров системы в 2 раза не приводит к качественному изменению кривых, представленных на рис. 5, а напряжение Пайерлса-Набарро оп-н изменяется не более чем на 1 МПа.
0.8 1 Время, нс
Рис. 5. Сдвиговое напряжение оху и соответствующее ему по времени положение дислокации ДD для модели после оптимизации с параметрами L = 10 нм, Н = 40 нм, АХ = 0,001 нм, Аt = 5 пс
ЗАКЛЮЧЕНИЕ
Полученные результаты молекулярно-динамического моделирования позволяют предположить, что для кристалла ОЦК-железа, ориентированного по направлениям [111], [0-11], [-211], содержащего краевую дислокацию, размер системы должен составлять не менее L = 10 нм, H = 40 нм, что соответствует системе, содержащей ~ 50000 атомов. При этом, величина сдвига атомов области В за один цикл не должна превышать 0,001 нм, а время одного цикла должно быть не менее 5 пс. Вычислительные мощности современных персональных компьютеров позволяют легко реализовать такие расчеты.
Величина напряжения Пайерлса-Набарро, полученная в настоящей работе сП-Н = 75 МПа хорошо согласуется с расчетами Hideki Mori и др. [6], где с помощью метода теории функционала плотности (density functional theory) и метода погруженного атома (потенциалы [7]) получены величины 80 МПа и 98 МПа соответственно.
Работа выполнена при поддержке РФФИ (грант 13-01-96051).
СПИСОК ЛИТЕРАТУРЫ
1. Нагорных И.Л., Бурнышев И.Н. Молекулярно-динамическое исследование декогезионного воздействия водорода на а-железо // Химическая физика и мезоскопия. 2011. Т. 13, № 1. С. 82-87.
2. Нагорных И.Л., Бурнышев И.Н. Численное моделирование влияния водорода на поведение кристаллов Al, Fe, Ni и Pd // Химическая физика и мезоскопия. 2012. Т. 14, № 4. С. 82-87.
3. Liu X., Golubov S.I., Woo C.H., Huang H. Glide of edge dislocations in tungsten and molybdenum // Materials Science and Engineering. 2004. A365. P. 960-100.
4. Johnson R.A., Oh D.J. Analytic embedded atom method model for bcc metals // Journal of Materials Research. 1989. V. 4, № 5. P. 1195-1201.
5. Wen M., Xu X-J., Fukuyama S, Yokogawa K. Embedded-atom-method functions for the body-centered-cubic iron and hydrogen // Journal of Matererials Research. 2001. V. 16, № 2. P. 3496-3502.
6. Mori H, Kimizuka H., Ogata Sh. Dislocation properties and Peierls stress of bcc iron based on generalized stacking fault energy surface by using first principles calculations // J. Japan Inst. Metals. 2009. V. 73, № 8. P. 595-600.
7. Mendelev M.I., Han S., Srolovitz D.J., Ackland G.J., Sun D.Y., Asta M. Development of new interatomic potentials appropriate for crystalline and liquid iron // Philosophical Magazine. 2003. V. 83, № 35. P. 3977-3994.
ABOUT MOTION SIMULATION EDGE DISLOCATIONS IN BCC-IRON
Nagornykh I.L., Burnyshev I.N.
Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. In the present work the optimization algorithm for BCC-Fe crystal model containing edge dislocation
a < 111 > {110} with periodical boundary conditions and EAM interatomic interactions is presented. Such optimization
allows to choose the model parameters at which the interaction between dislocation and its periodical image is negligible and the number of atoms in the modeling systems does not exceed 50000. The estimated value of Peierls-Nabarro stress is 75 MPa.
KEYWORDS: molecular dynamics, simulation, the edge dislocation, the Burgers vector, the Peierls-Nabarro stress.
Нагорных Иван Леонидович, кандидат физико-математических наук, научный сотрудник ИМ УрО РАН
Бурнышев Иван Николаевич, кандидат технических наук, ведущий научный сотрудник ИМ УрО РАН, тел. (3412) 20-74-33, e-mail: [email protected]
