Физика металлов и металловедение, 2021, T. 122, № 11, стр. 1164-1170

ВВЕДЕНИЕ

Молекулярно-динамические (МД) данные о первичной радиационной повреждаемости в ОЦК-металлах Fe и V свидетельствуют, что около 20% собственных межузельных атомов (СМА), образующихся в каскадах атомных столкновений с повреждающей энергией от 1 до 50 кэВ, объединены в димежузлия [1–4]. В матрице кристалла радиационные дефекты (РД) создают вокруг себя поле упругих деформаций, характеризуемое тензором дипольных моментов P_ij [5, 6]. Компоненты P_ij для димежузлий примерно в два раза больше, чем для одиночных СМА [7]. Это приводит к более сильному упругому взаимодействию димежузлий с внешними упругими полями по сравнению с одиночными СМА, т.к. в линейном приближении энергия взаимодействия равна [5]:

где S_ijkl – тензор упругих податливостей кристалла, а ε_ij – тензор упругих деформаций, создаваемых прилагаемым к кристаллу внешним полем напряжений σ_ij (здесь и далее используется правило Эйнштейна о суммировании по повторяющимся индексам). Таким образом, вклад димежузлий в радиационное деформирование кристалла под каскадообразующим облучением может быть заметным несмотря на то, что их образуется меньше, чем одиночных СМА. Оценка вкладов разных механизмов в деформирование кристалла при радиационном воздействии под нагружением (радиационные ползучесть, рост) показало, что наибольшим по величине эффектом обладает анизотропия диффузии РД, индуцированная приложенным напряжением [8]. Благодаря этому эффекту возникает зависимость эффективностей дислокационных стоков (ДС) для РД от приложенных внешних нагрузок, что приводит к перераспределению баланса потоков СМА (как одиночных, так и объединенных в кластеры) и вакансий к различно ориентированным относительно внешней нагрузки дислокациям и, как следствие, к переползанию дислокаций и, соответственно, к деформированию кристалла. Эффективность ДС определяется как обратный квадрат диффузионной длины РД за время их жизни – от образования до гибели на стоке, нормированный на плотность дислокаций. ДС, у которого отсутствует дальнодействующее упругое взаимодействие с РД, будем называть невзаимодействующим линейным стоком (НЛС). Относительные изменения эффективностей ДС и НЛС под действием приложенной нагрузки совпадают при относительно небольших величинах нагрузки (до 200–300 МПа в зависимости от температуры), как было показано в [9] кинетическим методом Монте-Карло (КМК). Это позволяет для описания влияния внешних нагрузок на эффективность ДС использовать полученное для НЛС аналитическое решение Саралидзе [8, 10]: где ξ₀(σ_ij) – эффективность НЛС под действием внешнего поля напряжений σ_ij, ξ₀(0) – эффективность НЛС в отсутствие внешнего поля напряжений, D_ij – тензор диффузии РД, TrD_ij – след тензора D_ij, D_⊥ – средний коэффициент диффузии в плоскости, перпендикулярной направлению НЛС. Таким образом, для определения изменения эффективностей ДС необходимо знать зависимости D_ij от внешнего упругого поля σ_ij или ε_ij.

D_ij при малых деформациях можно записать как степенное разложение по деформациям, ограничившись линейным членом [11]:

где D₀ – треть следа тензора D_ij в отсутствие внешних деформаций, d_ijkl – упругодиффузионный тензор, δ_ij – символ Кронекера. Тензор d_ijkl имеет ту же симметрию, что и тензор упругих постоянных кристалла c_ijkl, а значит, имеет три независимых константы d₁₁, d₁₂, d₄₄ (в нотации Фогта) для кристаллов кубической сингонии. В качестве трех независимых констант удобней выбрать три собственных значения тензора d_ijkl: d⁽¹⁾ = d₁₁ + 2d₁₂, d⁽²⁾ = d₁₁ – d₁₂, d⁽⁴⁾ = 2d₄₄. Совместное использование решения Саралидзе (2) и линейного разложения тензора D_ij по напряжениям или деформациям (в линейном приближении нет разницы каким именно разложением пользоваться, так как напряжения и деформации связаны линейно законом Гука) описывает эффективность НЛС с точностью не менее 1% относительно более точных значений КМК-метода, не использующего разложение (3), при |σ/T| ≤ 0.34 МПа/K в Fe и V (здесь T – температура, а σ – максимальная компонента тензора напряжений) [9]. Таким образом, для температурного интервала 300–1000 K разложение (3) можно считать выполняющимся с высокой точностью для деформаций меньших (1–4) × 10^–3 (разброс обусловлен выбором материала и конкретного значения температуры).

Тензор деформаций можно разложить разными способами. Для наших целей наиболее удобным является разложение [11] (здесь и далее используется кристаллографическая система координат с осями вдоль 〈100〉, если не оговаривается иное):

где

Соотношение (3) с учетом (4), (5) запишется как

Таким образом, для полного определения d_ijkl достаточно рассмотреть всего три определенных схемы нагружения: например, V¹, V², V⁴. Зная d_ijkl, можно определить значения дипольного тензора седловых конфигураций дефекта [11], позволяющие моделировать диффузию дефектов в неоднородных упругих полях.

Экспериментальное определение зависимости D_ij от внешнего упругого поля практически невозможно из-за значительных погрешностей [11], поэтому для определения таких зависимостей можно воспользоваться МД-методом.

Диффузионные характеристики димежузлий в Fe и V в отсутствие внешних упругих полей были получены нами МД-методом [12, 13]. В настоящей работе зависимости диффузионных характеристик димежузлий от величины однородной упругой деформации кристалла (в интервале от ‒1 до +1%) для деформированных состояний V¹, V², V⁴ рассчитывали МД-методом для ОЦК-металлов Fe (350–800 К) и V (350–600 К). Верхняя граница рассмотренного диапазона температур в Fe выше, чем в V, так как димежузлия в Fe более устойчивы к термической диссоциации по сравнению с V и, по большей части, успевают достичь микроструктурных стоков (дислокационная сетка плотностью 3 × 10¹⁴ м^–2) до диссоциации при температурах ниже 700–800 К [12]. По полученным деформационным зависимостям тензора диффузии димежузлий определяли упругодиффузионные тензоры димежузлий и эффективные дипольные тензоры седловых конфигураций димежузлий в Fe и V.

МЕТОДИКА МОДЕЛИРОВАНИЯ

Молекулярно-динамическая модель. МД-расчеты проводили с помощью программы для ЭВМ, разработанной авторами настоящей работы. МД-методики моделирования диффузионных траекторий дефектов в настоящей работе близко следуют методикам, описанным в обзорной работе [14]. Моделирование проводили с использованием потенциалов межатомного взаимодействия [7] и [15] для Fe и V соответственно, поскольку они хорошо описывают набор экспериментальных данных об объемных свойствах кристаллов и их собственных точечных дефектах [16]. Потенциалы разработаны в рамках подхода, в котором полная энергия модельного кристаллита представляется в виде суммы члена, зависящего от объема кристаллита, и члена, определяемого парными межатомными взаимодействиями. Модельный кристаллит представлял собой микроканонический ансамбль (постоянные число частиц, объем и полная энергия) с периодическими граничными условиями. Постоянную кристаллической решетки a до введения деформаций выбирали таким образом, чтобы при заданной температуре Т давление P в кристаллите было равно нулю с точностью до 0.1 эВ/нм³. Величину временного шага выбирали таким образом, чтобы среднее перемещение атомов за одну итерацию составило ~0.005a. Интегрирование уравнений движения выполняли методом Верле [17]. Положения двух СМА, составляющих димежузлие, определяли путем анализа заселенности ячеек Вигнера–Зейтца (ЯВЗ) атомами кристаллита. В качестве положения димежузлия использовалось положение одного из этих двух СМА.

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

При моделировании деформаций типа V¹ и V² использовали кристаллиты, содержащие N = 4394 + + 2 атомов, с поверхностями типа {100}. При моделировании деформаций типа V⁴ использовали кристаллиты, содержащие N = 4212 + 2 атомов, с поверхностями (110), $(1\bar {1}0),$ (001). При моделировании деформаций типа V¹, длины ребер кристаллита задавали равными 13(1 + ε₁)a. При моделировании деформаций типа V², длины ребер вдоль направлений [100], [010], [001] задавали равными 13(1 + ε₂)a, 13(1 – ε₂)a и 13a соответственно. При моделировании деформаций типа V⁴, длины ребер вдоль направлений [110], $[1\bar {1}0]$ и [001] задавали равными $9\sqrt 2 \left( {1 + {{\varepsilon }_{4}}} \right)a,$ $9\sqrt 2 \left( {1 - {{\varepsilon }_{4}}} \right)a$ и 13a соответственно. Величины ε₁, ε₂, ε₄ задавали равными 0, ±0.001, ±0005, ±0.01.

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

где x_j – координаты атома до деформации.

Для каждой температуры была произведена оценка времени, которое потребуется димежузлию, чтобы совершить 10⁵ скачков, на основании ранее полученных данных для среднего физического времени моделирования без деформации [13]. Это время служило критерием остановки моделирования каждой задачи.

Для каждого деформированного состояния при заданной температуре моделировали 100 диффузионных траекторий. Среднее число скачков, смоделированных для заданных температур и деформаций, составляло ∼6 × 10⁶, а среднее физическое время моделирования – от 2.1 до 4.3 мкс. Такой объем полученных данных позволил вычислять диффузионные характеристики с настолько высокой точностью, что влияние на них деформаций кристалла становится заметно даже при уровне деформаций в 0.1%.

Рассчитываемые характеристики. Для расчёта тензора диффузии дефектов D_ij каждую смоделированную диффузионную траекторию дефекта разбивали на некоторое число изохронных интервалов длительностью τ. Величину τ при каждой температуре и деформации выбирали исходя из компромисса между необходимостью в большом числе изохронных интервалов для точности расчета и необходимостью того, чтобы внутри изохронных интервалов содержалось достаточное число скачков для того, чтобы воспроизвести основные корреляции и физические особенности диффузии дефекта [14]. Для каждого изохронного интервала определяли вектор перемещения димежузлия R_nm, где n – номер траектории, а m – номер изохронного интервала. Затем вычисляли величину D = 〈R_nm ⊗ R_nm〉/2τ, где D – матрица, определяющая тензор D_ij в кристаллографических координатах, а операция 〈…〉 означает усреднение по всем n и m. Коэффициент диффузии D определяли как TrD_ij/3. Погрешности D и компонент тензора D_ij рассчитывали, как стандартные ошибки по указанному набору n и m.

Стартовая конфигурация димежузлия оказывает заметное влияние на рассчитываемый тензор диффузии димежузлия в Fe в отличие от V. Для того чтобы избавиться от этого влияния, при расчете диффузионных характеристик выбрасывали из обработки начальную часть каждой траектории, на которой имеется корреляция направлений миграции димежузлия с его стартовой ориентацией. Такая процедура эффективна для T ≥ 500 К. При более низких температурах эта корреляция сохранялась в течении всего времени моделирования траекторий, поэтому корректно определить тензоры диффузии димежузлий в Fe для T < 500 K не удалось, и значения их не приводятся.

Тензор диффузии для частных случаев. В процессе диффузии димежузлия в ОЦК-кристаллах Fe и V составляющие его СМА практически всегда совершают скачки в 1-е ближайшие соседи (БС). Доля скачков во 2-е БС составляет сотые доли процента от общего числа при всех рассматриваемых температурах и деформациях. Поэтому далее полагается, что димежузлие совершает скачки только в 1-е БС. Если димежузлие совершает случайные блуждания с вероятностью n_k совершить скачок в направлении k (k = 1, 2, 3, 4 для направлений скачков, параллельных $[111]$, $[\bar {1}11],$ $[1\bar {1}1],$ $[11\bar {1}]$ соответственно, n₁ + n₂ + n₃ + n₄ = 1), нормированный тензор диффузии ${{\tilde {D}}_{{ij}}}$ имеет вид

где A_ij = 0 при i = j и A_ij = n_i + 1 + n_j + 1 при i ≠ j. В табл. 1 приведены следующие из (8) частные случаи вида тензора ${{\tilde {D}}_{{ij}}}$ для конкретных значений n_k, которые пригодятся ниже.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Деформации типа V¹. На рис. 1 приведены зависимости D(ε₁) в Fe и V. Здесь и далее погрешности на рисунке не указываются, если они меньше размера символов на графике. Зависимость D от ε₁ заметно сильнее в V, чем в Fe, что свидетельствует о малости объема миграции димежузлия в Fe по сравнению с V. По всей видимости, это обусловлено тем, что в Fe значительный вклад в диффузию димежузлия вносит краудионный механизм диффузии в отличие от V [13]. Для Fe и V при всех значениях T и ε₁ выполняется $\tilde {D}_{{ij}}^{{}} = {{\delta }_{{ij}}}$ в пределах расчетных погрешностей.

Деформации типа V². На рис. 2 представлены зависимости D(ε₂) при 500 и 800 К в Fe и V. Зависимости D(ε₂) – четные, так как вследствие симметрии кристалла его поворот на 90° вокруг оси [001] эквивалентен изменению знака деформированного состояния на противоположный, а поворот кристалла не может изменить свойства его дефектов. Для Fe и V при всех значениях T и ε₂ выполняется $\tilde {D}_{{ij}}^{{}} = {{\delta }_{{ij}}}$ в пределах расчетных погрешностей.

Деформации типа V⁴. На рис. 3 представлены зависимости D(ε₄) в Fe и V для различных температур. Эти зависимости являются четными вследствие симметрии кристалла, как и в случае деформации типа V². При всех рассмотренных значениях T и ε₄ в пределах погрешности выполняются соотношения ${{\tilde {D}}_{{11}}} = {{\tilde {D}}_{{22}}}$ = ${{\tilde {D}}_{{33}}} = 1,$ ${{\tilde {D}}_{{23}}} = {{\tilde {D}}_{{13}}} = 0.$ Компонента ${{\tilde {D}}_{{12}}}$ ведет себя как нечетная функция (обусловлено симметрией кристалла), меняясь в пределах от –1 до 1 (рис. 4). Тензор ${{\tilde {D}}_{{ij}}}$ имеет вид ${{\tilde {D}}_{{ij}}} = {{\delta }_{{ij}}}$ ± ${{\delta }_{{i1}}}{{\delta }_{{j2}}} \pm {{\delta }_{{i2}}}{{\delta }_{{j1}}}$ при |ε₄| ≥ 1% (“+” для ε₄ ≥ ≥ 1%, “–” для ε₄ ≤ –1%), что, как видно из табл. 1, соответствует двумерной диффузии димежузлий в плоскостях $(\bar {1}10)$ и (110) при ε₄  ≥ 1% и ε₄  ≤ –1% соответственно.

Упругодиффузионные и дипольные тензоры. Собственные значения d^(α) (α = 1, 2, 4) для различных T получены путем аппроксимации МД-зависимостей D_ij (ε_α) соотношениями (6) в области малых деформаций (ε_α ∈ [–10^–3; 10^–3]). Полученные температурные зависимости d^(α) хорошо описываются выражениями: d⁽¹⁾ = –0.15D₀β, d⁽²⁾ = 0, d⁽⁴⁾ = 20.64D₀β в Fe; d⁽¹⁾ = –2.25D₀β, d⁽²⁾ = 0, d⁽⁴⁾ = = 12.80D₀β в V (β = (k_BT)^–1, где k_B – постоянная Больцмана).

В табл. 2 приведены полученные в [11] соотношения для кристаллов кубической сингонии, связывающие d^(α) с дипольными тензорами дефекта в основной $P_{{ij}}^{{\text{e}}}$ и седловой $P_{{ij}}^{{\text{s}}}$ конфигурациях при различной симметрии седловой конфигурации. В приведенных соотношениях используется обозначение $P_{0}^{{\text{e}}} = \frac{1}{3}{\text{Tr}}P_{{ij}}^{{\text{e}}}.$ Величины $P_{0}^{{\text{e}}}$ рассчитаны нами методом молекулярной статики: 36.11 эВ для Fe, 39.13 эВ для V.

Диффузия димежузлия происходит по разным механизмам, поэтому димежузлие обладает множеством различных седловых конфигураций. С помощью соотношений табл. 2 можно получить дипольный тензор “эффективной” седловой конфигурации, в который вносят вклад все реальные седловые конфигурации (каждая со своим весом). Из соотношения выражений в табл. 2 и полученных значений d^(α) видно, что для ОЦК-кристаллов Fe и V симметрия эффективной седловой конфигурации димежузлия является тригональной. Тогда получим: $P_{{ij}}^{{\text{s}}}$ = 36.06 эВ при i = j и $P_{{ij}}^{{\text{s}}}$ = = 10.32 эВ при i ≠ j для Fe; $P_{{ij}}^{{\text{s}}}$ = 38.38 эВ при i = j и $P_{{ij}}^{{\text{s}}}$ = 6.40 эВ при i ≠ j для V.

Знание дипольного тензора седловой конфигурации дефекта позволяет моделировать КМК-методом анизотропную диффузию дефекта в механических полях, создаваемых внешними нагружениями кристалла и/или внутренними источниками (дислокации, субграницы, др.), учитывая как эффекты упругой анизотропии кристаллов, так и симметрии кристалла [18]. Используемые в КМК-методе вероятности n_k могут быть представлены в виде [18]:

где $E_{k}^{{\operatorname{int} }}$ рассчитываются с помощью (1).

Продемонстрируем точность такого подхода на примере диффузии димежузлий в поле упругих деформаций типа V⁴. Подставив в уравнение (8) для компоненты ${{\tilde {D}}_{{12}}}$ выражения (1) и (9), получим

На рис. 5 представлены зависимости ${{\tilde {D}}_{{12}}}\left( {{{\varepsilon }_{4}}\beta } \right),$ рассчитанные МД-методом, а также с помощью выражений (6) и (10). Видно, что выражение (10), использующее в качестве параметра дипольный тензор, с высокой точностью описывает МД-данные для всех рассмотренных значений температур и деформаций, что выгодно отличает его от выражения (6), использующего в качестве параметра упругодиффузионный тензор, которое с хорошей точностью описывает МД-данные только при |ε₄β| < 0.05 эВ^–1.

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

ЗАКЛЮЧЕНИЕ

1. Разработаны методики МД-моделирования диффузии радиационных дефектов в ОЦК-кристаллах, подверженных нагрузкам различного типа (гидростатические, сдвиговые), и определения тензора диффузии по рассчитанным диффузионным траекториям. Эти методики использованы для расчета влияния внешних гидростатических и сдвиговых деформаций (в диапазоне значений от –1 до +1%) на тензор диффузии димежузлий в ОЦК-металлах Fe и V в температурном диапазоне 350–800 К.

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

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

Работа выполнена при частичной финансовой поддержке РФФИ (научный проект 18-08-01205-а) и НИЦ “Курчатовский институт” (приказ № 1934a от 28.09.2020) с использованием оборудования Центра коллективного пользования “Комплекс моделирования и обработки данных исследовательских установок мега-класса” НИЦ “Курчатовский институт” (http://ckp.nrcki.ru/).