Физика металлов и металловедение, 2020, T. 121, № 1, стр. 18-24

ВВЕДЕНИЕ

Высокая температура эксплуатации разрабатываемых энергетических реакторов четвертого поколения по сравнению с действующими водо-водяными энергетическими реакторами второго и третьего поколения [1] позволяет практически вдвое (до ≈ 60%) увеличить эффективность конверсии тепловой энергии в механическую за счет комбинирования термодинамических циклов Брайтона и Ранкина [2]. Однако повышение рабочей температуры реакторов зависит от создания новых конструкционных материалов, способных функционировать в условиях одновременного интенсивного термического воздействия, приложенных напряжений, в химически активной/окислительной среде под действием реакторного облучения. Благодаря своей исключительной по сравнению с аустенитными и мартенситно-ферритными сталями коррозионной стойкости и сопротивлению термической ползучести, жаропрочные никелевые сплавы рассматриваются в качестве основных конструкционных материалов для использования в трех из шести существующих дизайнов ядерных энергетических установок четвертого поколения [1, 3].

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

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

ФОРМУЛИРОВКА ЗАДАЧИ

Каскады смещений инициируются ПВА с энергиями Е_ПВА ≥ 1 кэВ и являются основным источником радиационных повреждений, создаваемых в материалах, подвергаемых облучению быстрыми частицами в режиме упругих потерь энергии. Оценка характерных линейных размеров и времени релаксации каскадов смещений составляет ≈5‒30 нм и ≈2–20 пс соответственно в зависимости от материала мишени, энергии ПВА и температуры облучения. Малые времена и размеры не позволяют исследовать процессы, протекающие в каскадах, экспериментальными методами, но первичное дефектообразование в каскадах смещений может быть смоделировано методом молекулярной динамики (МД).

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

ИСПОЛЬЗОВАННЫЕ МЕТОДЫ И ПОДХОДЫ

Для вычисления сил межатомного взаимодействия в никеле использован межатомный потенциал [8], построенный по методу внедренного атома. На коротких расстояниях парная часть потенциала модифицирована подстановкой универсального потенциала Зиглера–Бирсака–Литтмарка [9]. Экспериментально измеренная пороговая энергия смещения в никеле E_d = 23 ± 2 эВ [10, 11] использована в качестве подгоночного параметра. Пороговая энергия смещения модифицированного потенциала составляет 23 < E_d ≤ 24 эВ. Для подгонки потенциала использовали процедуру, описанную в [12]. Подгонка потенциала никак не повлияла на равновесный параметр решетки, энергию когезии E₀ = –4.45 эВ, энергию образования вакансии $E_{{\text{y}}}^{{\text{f}}}$ = 1.57 эВ, модули упругости, энергию дефектов упаковки, энергию свободной поверхности и т.п.

Каскады смещений в никеле смоделированы для температур T = 100, 300, 600, 900 и 1200 K. Теорема о вириале [13] использована для определения температурной зависимости равновесных параметров решетки, соответствующих нулевым внутренним напряжениям. Значения равновесных параметров решетки при температурах моделирования приведены в табл. 1.

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

Каскады смещений инициированы ПВА с энергией Е_ПВА = 5, 10, 15 и 20 кэВ. Чтобы смоделировать изотропное пространственное и случайное временно́е распределение ПВА, их вводили в разных местах кристалла вдоль одного из кристаллографических направлений 〈123〉 в различные моменты времени. Для каждой пары значений параметров (Е_ПВА, T) была смоделирована серия из 24 каскадов смещений.

Размер моделируемого кристалла масштабировали в зависимости от энергии ПВА в пропорции ≈10^–2 эВ/атом (см. табл. 2) так, чтобы каскады смещений не пересекали границы кристалла. Перед введением ПВА кристалл никеля приводили в состояние термодинамического равновесия для температуры моделирования в течение 1 × 10⁴ МД-итераций. МД-моделирование проводили без контроля температуры. Пример характерного изменения эффективной температуры Максвелла на разных стадиях эволюции каскада смещений, инициированного ПВА с энергией Е_ПВА = 20 кэВ, в никеле при температуре облучения T = 100 K показан на рис. 1. Энергия, вносимая ПВА, не извлекалась из системы, а соответствующее повышение температуры кристалла после релаксации каскада смещений не превышало ≈40 град ни в одном из смоделированных кристаллов.

Рис. 1.

Эффективная температура Максвелла T, шаг интегрирования по времени τ и число пар Френкеля N_FP на различных этапах эволюции каскада смещений, инициированного ПВА с энергией Е_ПВА = 20 кэВ в никеле при температуре 100 K.

На начальной стадии развития каскада смещений относительно небольшое число атомов кристалла движется с большой скоростью, в то время как основной объём материала продолжает находиться в состоянии термодинамического равновесия. В методе скоростей Верле [14], использованном для интегрирования уравнений движения в этой работе, для сходимости решения шаг интегрирования по времени τ выбирается, исходя из энергии самого быстрого атома. Сразу после введения ПВА шаг интегрирования τ, обеспечивающий сходимость, падает на три порядка величины (см. рис. 1). Таким образом, прямое интегрирование уравнений движения всего ансамбля ведет к неэффективному использованию вычислительных ресурсов. Для оптимизации вычислений на начальной стадии развития каскадов смещений использован метод [15]. Устойчивость алгоритма [15] протестирована ранее, а сам метод неоднократно применяли при моделировании радиационных повреждений в алюминии [16, 17], меди [18, 19], α-цирконии [19, 20], интерметаллидах γ‑TiAl [21], α₂-Ti₃Al [22] и γ'-Ni₃Al [23], а также при исследованиях взаимодействия каскадов смещений с дислокациями [24, 25].

Для идентификации и визуализации радиационных дефектов использованы критерий Линдеманна [26], метод ячеек Вигнера–Зейтса [27] и кластерный анализ [18]. Пороговый радиус, использованный в критерии Линдеманна, равен 0.3a, где a – это равновесный параметр решетки. Радиус первой координационной сферы a/√2 применяли при определении кластеров точечных дефектов в кластерном анализе.

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

Число пар Френкеля N_FP, образованных в индивидуальных каскадах смещений в никеле, и соответствующие средние значения 〈N_FP〉 как функция (Е_ПВА, T) приведены на рис. 2. Зависимость 〈N_FP〉 от Е_ПВА близка к линейной при всех смоделированных условиях облучения. Если энергия ПВА задана в единицах кэВ, среднее число пар Френкеля аппроксимируется соотношением 〈N_FP〉 = 2 ± 0.9 × $E_{{{\text{ПВА}}}}^{{1.1 \pm 0.1}}$. Простая функциональная форма позволяет легко экстраполировать зависимость 〈N_FP(Е_ПВА)〉 в область Е_ПВА > 20 кэВ. С увеличением Е_ПВА наиболее резкий рост 〈N_FP〉 происходит при низкой температуре облучения. Увеличение температуры ведет к постепенному снижению 〈N_FP〉, и наиболее резкое падение наблюдается при высоких энергиях ПВА.

Рис. 2.

Число пар Френкеля N_FP, образованных в каскадах смещений в никеле. Открытыми символами показаны значения N_FP после релаксации отдельных каскадов (показано со смещением). Соответствующие средние значения 〈N_FP〉 в серии каскадов с одинаковыми параметрами (Е_ПВА, T) показаны закрашенными символами. Вертикальными отрезками здесь и далее обозначены стандартные ошибки среднего.

Дисперсия значений N_FP велика даже для каскадов смещений, смоделированных при одинаковых значениях параметров (Е_ПВА, T). Разброс N_FP увеличивается с ростом Е_ПВА и снижением T. Чтобы получить статистически достоверные значения N_FP в зависимости от (Е_ПВА, T), следует определить минимально необходимый размер статистической выборки. В проводимом исследовании для определения размера выборки использована зависимость N_FP от числа n смоделированных каскадов в серии с одинаковыми значениями (Е_ПВА, T) (см. пример на рис. 3). С ростом n, 〈N_FP〉 сходится к своему “стационарному” значению, определяя минимально необходимый размер ñ выборки (в примере на рис. 3 соответствует соотношению ñ ≥ 18).

Рис. 3.

Зависимость 〈N_FP〉 от числа n каскадов в статистической выборке. Каскады инициированы ПВА с энергией Е_ПВА= 20 кэВ.

Доля σ_vac =Σ N_vac/N_FP вакансий в вакансионных кластерах размером N_vac ≥ 3 и σ_SIA =Σ N_SIA/N_FP междоузлий (SIAs) в междоузельных кластерах размером N_SIA ≥ 4 определены и усреднены по серии каскадов с одинаковыми значениями (Е_ПВА, T), рис. 4. Как 〈σ_vac〉, так и 〈σ_SIA〉 демонстрируют похожую зависимость от Е_ПВА при всех условиях моделирования. Единственное исключение – это постоянное значение 〈σ_vac〉 ≈ 0.06–0.07 при T = 1200 K, т.е. фактическое отсутствие кластерообразования и доминирование изолированных вакансий и дивакансий в облученном никеле при высоких температурах.

Рис. 4.

Доля вакансий 〈σ_vac〉 и междоузлий 〈σ_SIA〉 в кластерах точечных дефектов в зависимости от (Е_ПВА, T).

Зависимость 〈σ_vac〉 повторяет зависимость среднего числа вакансионных кластеров 〈Y_vac〉 на каскад (см. рис. 5 для сравнения). Число создаваемых радиационных дефектов в 5 кэВ каскадах смещений мало́ при всех температурах облучения за исключением T = 100 K (см. рис. 2). По этой причине как 〈σ_vac〉, так и 〈Y_vac〉 малы при 300 ≤ T ≤ 1200 K и относительно велики при низкой температуре.

Рис. 5.

Зависимость среднего числа вакансионных кластеров на каскад 〈Y_vac〉 от условий моделирования (Е_ПВА, T).

Увеличение 〈N_FP〉, вызванное увеличением Е_ПВА от 5 до 10 кэВ, сопровождается увеличением 〈Y_vac〉 и 〈σ_vac〉 при температуре 300 ≤ T ≤ 900 K. Такое же увеличение энергии ПВА при T = 100 K не приводит к росту 〈Y_vac〉 в силу того, что число создаваемых вакансий просто недостаточно для образования более одного кластера. Однако дальнейшее увеличение энергии ПВА и, следовательно, N_FP ведет к устойчивому росту как 〈Y_vac〉, так и 〈σ_vac〉 в каскадах смещений при T = 100 K.

Величины 〈Y_vac〉 и 〈σ_vac〉 принимают одни и те же значения при моделировании каскадов смещений, инициированных ПВА с энергиями 10  ≤ Е_ПВА ≤ 15 кэВ при температуре моделирования в диапазоне 300 ≤ T ≤ 900 K. Согласно рис. 2, число 〈N_FP〉, создаваемое ПВА с этими энергиями при температуре 300 ≤ T ≤ 900 K приблизительно равно числу 〈N_FP〉, создаваемому ПВА с энергиями 5 ≤ Е_ПВА ≤ 10 кэВ, но при температуре T = 100 K. Таким образом, в каскадах смещений в никеле зависимость 〈Y_vac〉 от 〈N_FP〉 имеет три стадии: рост с увеличением 〈N_FP〉 при низких значениях 〈N_FP〉, стационарное состояние при промежуточных значениях 〈N_FP〉, когда число образованных дефектов недостаточно для образования более одного вакансионного кластера, и дальнейший рост при больших значениях 〈N_FP〉. Доля вакансий в вакансионных кластерах 〈σ_vac〉 повторяет функциональную зависимость 〈Y_vac〉 от (Е_ПВА, T).

Средний размер вакансионных кластеров 〈N_vac〉, зародившихся в каскадах смещений в никеле, как функция (Е_ПВА, T) показан на рис. 6. Три характерных типа зависимости 〈N_vac(Е_ПВА)〉 определяются температурой облучения и термической стабильностью вакансионных кластеров. При температуре T = 100 и 300 K вакансионные кластеры стабильны, 〈N_vac〉 ∝ 〈N_FP〉 и, следовательно, ∝ Е_ПВА. При температуре T = 1200 K вакансионные кластеры нестабильны, средний размер 〈N_vac〉 мал и практически не зависит от Е_ПВА. При температурах T = 600 и 900 K стабильность вакансионных кластеров определяется их размером и степенью совершенства.

Рис. 6.

Зависимость среднего размера вакансионных кластеров 〈N_vac〉 от параметров (Е_ПВА, T).

В никеле и никелевых сплавах экспериментально наблюдается два типа вакансионных кластеров (см, напр., [28, 29]). Маленькие вакансионные кластеры преимущественно образуют тетраэдры дефекта упаковки, которые с ростом вакансионных кластеров превращаются в вакансионные петли Франка с вектором Бюргерса 1/3〈111〉. Термическая стабильность тетраэдров дефекта упаковки сильно зависит от их размера и уровня совершенства. Тетраэдры дефекта упаковки правильной формы с числом вакансий N_vac = (k + 1)k/2, k = 2, 3, 4… непропорционально более стабильны по сравнению с нерегулярными тетраэдрами дефекта упаковки, число вакансий в которых отклоняется от этих “магических чисел”. По этой причине при температуре T = 600 K в каскадах смещений, инициируемых ПВА с энергиями 5 ≤ Е_ПВА ≤ 15 кэВ и Е_ПВА = 20 кэВ, преимущественно зарождаются вакансионные кластеры размером 〈N_vac〉 = 6 и 10 соответственно, рис. 6. По той же самой причине вакансионные кластеры с 〈N_vac〉 ≈ 6 и 10 образуются в каскадах смещений при T = 900 K, однако их среднее число 〈Y_vac〉 оказывается вдвое меньше, чем в каскадах смещений при температуре T = 600 K (см. рис. 5).

Зависимость 〈σ_SIA〉 на рис. 4 совпадает с 〈σ_vac〉, однако в отличие от 〈σ_vac〉, 〈σ_SIA〉 определяется средним числом междоузлий 〈N_SIA〉 в междоузельных кластерах (рис. 7). При низких энергиях ПВА 〈N_SIA〉 ≈ 7. Это число соответствует размеру икосаэдрического междоузельного кластера в ГЦК-структуре [30]. Увеличение Е_ПВА ведет к увеличению N_FP и вызывает рост 〈N_SIA〉. Немонотонный рост 〈N_SIA〉 связан с образованием нескольких стабильных междоузельных кластеров размером 〈N_SIA〉 = 7, 10, 13, 15 и 19 междоузлий.

Рис. 7.

Зависимость среднего размера междоузельных кластеров 〈N_SIA〉 от (Е_ПВА, T).

Термоактивируемое превращение неподвижных междоузлий, имеющих конфигурацию гантели вдоль 〈100〉, в подвижную конфигурацию вдоль 〈110〉 подавлено при температуре T = 100 K. Из-за низких температур не наблюдается изменение размеров междоузельных кластеров в сторону ближайшей стабильной конфигурации, и 〈N_SIA〉 изменяется пропорционально как 〈N_FP〉, так и Е_ПВА.

Среднее число междоузельных кластеров на каскад 〈Y_SIA〉 в зависимости от (Е_ПВА, T) показано на рис. 8. Чтобы сделать корреляцию между 〈Y_SIA〉 и 〈N_FP〉 еще более очевидной, результаты моделирования, показанные на рис. 2, были перегруппированы и добавлены к зависимости 〈Y_SIA〉 на рис. 8, из которого видно, что 〈Y_SIA〉 ∝ 〈N_FP〉 при всех смоделированных условиях облучения (Е_ПВА, T).

Рис. 8.

Величины 〈Y_SIA〉 и 〈N_FP〉 как функции (Е_ПВА, T).

ВЫВОДЫ

Компьютерное моделирование методом МД использовано для исследования процесса первичного дефектообразования в каскадах смещений в никеле, подвергаемом облучению быстрыми частицами в режиме упругих потерь энергии. Были смоделированы каскады смещений, инициированные ПВА с энергиями 5 ≤ Е_ПВА ≤ 20 кэВ в материале при температуре 100 ≤ T ≤ 1200 K. Число пар Френкеля, доля вакансий и междоузлий в кластерах точечных дефектов, средний размер вакансионных и междоузельных кластеров, и среднее число кластеров на каскад определены как функция (Е_ПВА, T). Чтобы получить статистически достоверные количественные результаты, для каждой пары (Е_ПВА, T), смоделирована серия из 24 каскадов смещений.

Зависимость 〈N_FP〉 от энергии ПВА аппроксимируется близкой к линейной показательной функцией 〈N_FP〉 = 2 ± 0.9 × $E_{{{\text{ПВА}}}}^{{1.1 \pm 0.1}}$ при всех смоделированных параметрах (Е_ПВА, T). Установлено, что как 〈σ_vac〉, так и 〈σ_SIA〉 демонстрируют похожую функциональную зависимость от Е_ПВА, однако физические механизмы, определяющие их поведение, различны. Значение 〈σ_vac〉 повторяет функциональную зависимость 〈Y_vac〉 при всех условиях облучения (Е_ПВА, T). Величина 〈σ_SIA〉 определяется размером междоузельных кластеров 〈N_SIA〉 и диффузионной подвижностью междоузлий. При температуре T = 100 K мобильность междоузлий подавлена, и 〈σ_SIA〉 ∝ Е_ПВА. При T ≥ 300 K, междоузлия образуют как небольшие междоузельные дислокационные петли с вектором Бюргерса 1/2〈110〉, так и икосаэдрические кластеры с N_SIA = 7 [30] или кластеры, производные от икосаэдрических [31].

Размеры вакансионных кластеров 〈N_vac〉 определяются температурой облучения и термической стабильностью тетраэдров дефекта упаковки. При температурах T ≤ 300 K все вакансионные кластеры стабильны, 〈N_vac〉 ∝ Е_ПВА. При температурах 600 ≤ T ≤ 900 K стабильны преимущественно вакансионные кластеры с N_vac = 6 и 10, размер которых соответствует регулярным тетраэдрам дефекта упаковки. Значения 〈Y_SIA〉 ∝ 〈N_FP〉 и, следовательно, ∝ Е_ПВА при всех условиях облучения.

Исследования выполнены при поддержке НИЦ “Курчатовский институт”, проект № 1603. Программное обеспечение для моделирования радиационных эффектов методом МД, численные методы интегрирования и методы идентификации и визуализации дефектной структуры материалов разработаны при поддержке Российского фонда фундаментальных исследований, проект РФФИ № 17-03-01222a. Моделирование первичных повреждений выполнено с использованием высокопроизводительных вычислительных ресурсов центра коллективного пользования “Комплекс моделирования и обработки данных исследовательских установок мега-класса” НИЦ “Курчатовский институт”, http://ckp.nrcki.ru.