CC BY
116
УДК 621.039.531.001.57
М. Ю. Тихончев, В. В. Светухин, Т. С. Ильина
МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ПЕРВИЧНОЙ РАДИАЦИОННОЙ ПОВРЕЖДАЕМОСТИ а-ЖЕЛЕЗА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ*
В работе представлены результаты компьютерного моделирования процессов первичной радиационной повреждаемости а-железа методом молекулярной динамики с использованием расчетного кода FRENKLOW. При моделировании использован многотельный потенциал межатомного взаимодействия. Рассчитаны пороговые энергии смещения для различных кристаллографических направлений, промоделированы каскады смещений для энергий первично выбитого атома от 0,1 до 20 КэВ и проведены оценки количества дефектов, переживающих рекомбинацию в каскаде, получены результаты по количеству и размерам кластеров вакансий и внедрений, образующихся в таких каскадах.
Введение
Одной из важных характеристик радиационного воздействия, используемых в настоящее время при исследовании изменений физических и механических свойств материалов при длительном облучении, является так называемая повреждающая доза, которая выражается в количестве атомных смещений в кристаллической решетке материала на атом (сна). В физике радиационной повреждаемости материалов предложено немало различных математических моделей, для расчета повреждающей дозы, среди которых наибольшее признание и распространение получила модель Норжетта-Робинсона-Торренса (NRT-стандарт) [1]. Однако эта модель имеет ряд существенных недостатков, основным из которых является неоднозначность взаимосвязи количества смещений на атом с реальными изменениями микроструктуры материала при различных условиях облучения. В последние годы разрабатываются новые концепции, основанные на учете процессов рекомбинации точечных дефектов в каскадах и фракций свободно-мигрирующих дефектов (freely migrating defects, - FMD) [2, 3]. Они призваны более корректно оценивать повреждающую способность различных энергетических групп нейтронного спектра. Для развития таких концепций необходимо получение детальной информации об особенностях процессов первичного радиационного повреждения материалов, происходящих в каскадах атомных смещений. Вследствие скоротечности таких каскадов их исследование экспериментальными методами затруднено. Поэтому основным инструментом их изучения является математическое моделирование.
Данная работа посвящена моделированию каскадов атомных смещений в а-железе методом молекулярной динамики с целью определения ряда параметров первичной радиационной повреждаемости с учетом процессов рекомбинации и кластеризации точечных дефектов в каскадах смещений.
1. Молекулярно-динамическая модель
Для моделирования каскадов атомных смещений в а-железе нами использовалась компьютерная программа FRENKLOW [4], разработанная в
* Работа выполнена при поддержке гранта РФФИ 07-02-92282-СИГ_а РНП 2.1.2.7242. 70
Тверском государственном техническом университете. Программа базируется на методе молекулярной динамики и предназначена:
- для расчета динамическим методом пороговых энергий смещения и замещения для различных направлений вылета первично-выбитого атома (ПВА);
- для моделирования процессов развития цепочек атом-атомных соударений;
- для моделирования каскадов атомных смещений;
- для расчета атомных конфигураций точечных дефектов и их кластеров;
- для расчета энергий образования вакансий и миграции вакансий и межузельных атомов.
В настоящей работе использовался полуэмпирический многотельный потенциал межатомного взаимодействия. При использовании многотельных потенциалов энергия каждого атома не имеет вида суммы парных взаимодействий, а представляется как некоторая функция его локального окружения. Для переходных металлов были разработаны несколько схем построения такого потенциала: метод погруженного атома (МПА) [5], схема Финниса-Синкляра [6] и схема Росато-Гвиллопа-Легранда [7]. Несмотря на несколько различные физические интерпретации, все эти методы дают одинаковое аналитическое выражение для полной энергии системы из N частиц:
N -1 N N
Ео = х Е ф(гу)+X р(р«); (1)
1 =1 у=г+1 1=1
N
Рг = X ^(Гу ^ (2)
У=1 у *«
где Ео - полная энергия системы; в формализме МПА: ^р«) - функция внедрения, определяющая энергию внедрения атома в электронную жидкость плотностью р7; ^(гу) - собственная электронная плотность у-го атома как функция расстояния до его центра; ф(гу) - парный потенциал взаимодействия между атомами г и у; Гу - расстояние между атомами г и у.
В настоящем исследовании мы использовали функцию внедрения, функцию электронной плотности и парную часть потенциала для а-Бе из работы [8]. Проверка потенциала проведена путем расчета и сравнения с соответствующими экспериментальными значениями равновесного параметра решетки (а0), энергии связи (Ес), упругих констант (Сц, С\2, С44), энергии образования вакансии (Еу) с учетом и без учета релаксации кристаллита, энергии активации вакансии (Е0), энергии образования внедрения (Е1^ типа «гантель» для трех ориентаций: <100>, <110> и <111> (рис. 1). Полученные рассчитанные величины вместе с соответствующими экспериментальными значениями представлены в таблице 1.
Отметим, что параметры а0, Си, Сц, С44 и ЕУ]‘ использовались как подгоночные при построении потенциала. Поэтому по этим параметрам наблюдается очень хорошее согласие расчета и эксперимента. Для остальных рассмотренных параметров наибольшие отклонения от эксперимента наблюдаются для энергии активации вакансии (-11%), энергии образования «гантели» внедрения с ориентацией вдоль направления <100> (-10%) и энергии
связи (-7%). Отклонение для Е$ <110> и Е$ <111> составили +2,6 и +1,5% соответственно. В целом согласие между рассчитанными и соответствующими экспериментальными величинами следует признать удовлетворительным.
Рис. 1 Ориентации «гантелей» внедрений
Таблица 1
Экспериментальные и рассчитанные параметры а-железа
эксперимент расчет эксперимент расчет
ао, А Ес, эВ С11, ГПа С12, ГПа С44, ГПа Еу, эВ 2,8553 [9]; 2,866 [10] 4,316 [9]; 4,28 [14] 243,4 [9]; 233,0 [11] 145.0 [9]; 135,5 [11] 116.0 [9]; 117,8 [11] 1,84 [9]; 2,00 [12] 2,8553 4,013 243,3 145.0 116.1 1,84 1,72 (релакс.) Ев, эВ Е‘/ <100>, эВ Е‘/ <110>, эВ £%<111>,эВ Е27<110>, эВ Е2^ <100>, эВ 2,65 [15] 4,37 [13] 3,41 [13] 4,11 [13] 2,35 4,17 3,53 3,94 3.45 3,30 ( релакс .) 3.46 3,20 (релакс.)
В таблице 1 также представлены расчетные энергии образования бивакансии (Е1У) для двух ориентаций, а именно вдоль направлений <110> и <100>. Полученные значения энергии образования таких бивакансий практически совпадают, если не учитывать релаксацию кристаллита. При учете релаксации энергия для бивакансии вдоль направления <100> оказывается немного ниже (на = 3%), чем вдоль направления <110>.
На рисунке 2 представлены зависимости изменения энергии кристаллита от изменения ориентаций «гантелей» внедрения от направлений <111> и <100> к направлению <110>.
Рис. 2 Изменение энергии кристаллита при повороте «гантели» внедрения: а - от ориентации <111> к ориентации <110>; б - от ориентации <100> к ориентации <110>
При повороте от <111> к <110> (рис. 2,а) величина энергетического «барьера» составляет лишь =0,001 эВ, что находится в пределах погрешности расчета. Таким образом, ориентацию <111> при используемом потенциале межатомного взаимодействия можно считать неустойчивой. При повороте от <100> к <110> (рис. 2,б) величина аналогичного энергетического «барьера» составляет 0,045 эВ.
2. Пороговые энергии смещения
Для моделирования кристаллической структуры а -железа при расчетах пороговой энергии смещения (Ел) использовались ОЦК кристаллиты, содержащие около 10000 атомов. При расчетах использовались «периодические» граничные условия, которые позволяют имитировать бесконечную протяженность кристалла. При решении задачи нахождения пороговой энергии поиск устойчивых конфигураций начинается с небольших энергий ПВА, что приводит на начальном этапе расчета к образованию неустойчивых смещений. Пошагово повышая энергию ПВА, можно определить энергию, при которой в кристаллите образуются устойчивые пары Френкеля. При реализации этой процедуры мы повышали энергию с шагом 1 эВ. Таким образом, погрешность определения Ел для фиксированного кристаллографического направления, обусловленная дискретным изменением энергии ПВА, не превышает здесь 1 эВ. Следует заметить, что метод пошагового повышения энергии нельзя заменить методом бинарного поиска, поскольку наблюдение появления устойчивой пары Френкеля при некоторой энергии ПВА, вообще говоря, не гарантирует образования такой пары при большей энергии.
Согласно анализу, проведенному в работе [16], имеется небольшое количество экспериментальных результатов по определению пороговых энергий смещения для а-железа, полученных путем облучения электронами тонких монокристаллических фольг. Однако, согласно работам [17-19], в таких экспериментах трудно достичь надежных результатов. В частности низкоэнергетический порог для некоторых направлений может инициироваться электронами, движущимися по наклонной к этим направлениям. Таким образом, надежные результаты удается получить только для минимального значения Ел, которое составляет для а-железа 16-20 эВ и соответствует кристаллографическим направлениям вблизи <100>.
Рассчитанные нами значения Ел в сравнении с известными экспериментальными результатами представлены в таблице 2. Следуя методу, предложенному в [16], мы определяли пороговую энергию для заданного кристаллографического направления как минимальное из рассчитанных значений Ел для всевозможных направлений, отклоняющихся от заданного на угол до 20°. Полученные нами значения Ел хорошо согласуются с большинством экспериментальных результатов. Заметное расхождение наблюдается только для направления <111>. Отметим, что рассчитанные нами значения Ел несколько отличаются от результатов, полученных авторами работы [16], которые, используя тот же потенциал, получили 17, 33 и 33 эВ для направлений <100>, <110> и <111> соответственно и 15 эВ - как минимальное значение.
При расчете повреждающих доз важную роль играет величина средней пороговой энергии смещения Ел . Для оценки средней пороговой энергии Ел проведена серия расчетов для различных направлений импульса ПВА, кото-
рые задавали путем моделирования изотропного случайного вектора. Всего было промоделировано 400 различных направлений. Все расчеты проводились при начальной температуре кристаллита 0 К. По результатам проведенных расчетов получена оценка Е^ = 38,9 ± 1,5 эВ. Эта оценка близка к используемому в настоящее время значению Е^ = 40 эВ, рекомендованному ЛБТМ стандартом для определения повреждающих доз в железе и сталях [23].
Таблица 2
Рассчитанные и экспериментальные значения пороговых энергий смещения, эВ
Минимальное Кристаллографическое направление
значение <100> <110> <111>
расчет 18 18 29 28
Эксперимент [20] 17 17 >30 20
[21] 20 20 30 -
[22] 16-18 - - -
3. Моделирование каскадов смещений
Для моделирования каскадов атомных смещений и оценки числа «выживших» дефектов задавались ОЦК кристаллиты железа, содержащие до =600000 атомов. При этих расчетах также использовались «периодические» граничные условия. Моделирование проводилось при нулевой начальной температуре кристаллита для десяти различных энергий ПВА (ЕПВА): 0,1; 0,5; 1; 2; 3; 4; 5; 10; 15 и 20 КэВ. Расчеты проводились с неравномерным шагом по времени, который выбирался так, чтобы он не превосходил 10-3 пс, и чтобы за один шаг по времени атом с максимальной кинетической энергией смещался не более чем на 0,03 А.
Моделируемое время развития каскада подбиралось так, чтобы обеспечить моделирование всего процесса образования и релаксации дефектов в каскаде вплоть до его затухания. Информация по размерам модельного кристаллита и моделируемого времени для рассматриваемых энергий ПВА собрана в таблице 3.
Таблица 3
Размеры кристаллита и моделируемое время
Энергия ПВА, КэВ Число атомов кристаллита Моделируемое время, пс
0,1 =100 тыс 5
0,5 =150 тыс. 10
1 =200 тыс. 10
2 =200 тыс. 15
3 =250 тыс. 20
4 =300 тыс. 30
5 =300 тыс. 40
10 =400 тыс. 50
15 =500 тыс. 60
20 =600 тыс. 60
Для каждой из рассматриваемых энергий моделировалось девятнадцать направлений импульса ПВА:
- 5 направлений, равномерно распределенных между <100> и <110>;
- 5 направлений, равномерно распределенных между <110> и <111>;
- 5 направлений, равномерно распределенных между <001> и <111>;
- 4 направления, равномерно распределенные между <110> и <101>.
После завершения моделирования каскада смещений проводился анализ кристаллита, подсчитывалось число дефектов, переживших рекомбинацию в каскаде, и определялось среднее число таких дефектов для каждой энергии ПВА. Подсчет дефектов в кристаллите осуществлялся следующим образом. Каждому узлу I идеальной кристаллической решетки ставится в соответствие ячейка Дирихле С, которая является выпуклым многогранником и определяется как множество всех точек пространства, расстояние от которых до узла I (с учетом периодических граничных условий) меньше или равно расстоянию до любого другого узла решетки. Отсутствие атомов в узле С, свидетельствует о вакансии в этом узле, попадание более одного атома в ячейку С{ трактуется как наличие внедрений вблизи узла I. Число дефектов определяется как общее количество ячеек Дирихле, не содержащих ни одного атома материала.
Доля «выживших» дефектов определялась по формуле
р(ЕПВА ) _
N ( Епва ) /( епва )
(3)
где ЩЕша) - рассчитанное среднее число «выживающих» дефектов, /^ЕПВА) = = 0,8 • Епва/(2 ) - количество атомных смещений по модели Кинчина-Пиза.
Величину р(ЕПВА) в современной литературе часто называют коэффициентом каскадной эффективности. Для пороговой энергии смещения в модели Кинчина-Пиза мы использовали значение Е^ = 40 эВ. Полученные значения И(Епва) и р(ЕПВА) представлены в рисунках 3 и 4 соответственно.
Энергия ПВА, КэВ
Рис. 3 Рассчитанное число «выживающих» дефектов
Энергия ПВА, КэВ
Рис. 4 Доля дефектов, переживших рекомбинацию в каскаде смещений
Из полученных результатов видно, что везде на рассматриваемом интервале энергий, кроме диапазона от 3 до 5 КэВ, наблюдается рост числа де-
фектов с увеличением энергии. В диапазоне 3-5 КэВ наблюдается уменьшение числа дефектов с ростом энергии ПВА. Доля дефектов, переживших рекомбинацию, уменьшается с увеличением энергии ПВА до 10 КэВ, и при энергии 10 КэВ достигает значения 0,27. В диапазоне ЕПВА от 0,5 до 3 КэВ изменение р(ЕПВА) идет по закону, близкому к линейному. В интервале ЕПВА = = 3 - 5 КэВ наблюдается резкое снижение числа «выживающих» дефектов, которое продолжается до ЕПВА = 5 КэВ. При энергиях выше 10 КэВ функция р(ЕПВА) становится возрастающей и достигает значения 0,47. Такое поведение функции р(ЕПВА) можно объяснить тем, что при энергиях выше 10 КэВ начинает наблюдаться расщепление каскада на субкаскады. Энергия субкаскадов ниже энергии исходного каскада, и доля дефектов, переживающих рекомбинацию в них, оказывается выше.
Известно, что существенный вклад в микроструктурную эволюцию материала под облучением вносит объединение производимых в нем точечных дефектов в кластеры. При моделировании каскадов смещений мы, наряду с оценкой числа «выживающих» дефектов, получили оценки размеров и количества кластеров внедрений и вакансий, остающихся в кристаллите после затухания каскада. Дефекты одного типа считали принадлежащими одному кластеру, если границы соответствующих им ячеек Дирихле имели общий участок. Полученные результаты представлены на рисунке 5 и в виде гистограмм на рисунке 6. Хорошо видно, что при энергиях ПВА до 15 КэВ количество вакансий, участвующих в процессе кластеризации, заметно больше числа внедрений, входящих в кластеры. При ¿Пад. = 20 КэВ число вакансий в кластерах лишь немного превышает число внедрений, образующих кластеры. Начиная с £’ПВА = 5 КэВ вероятность образования кластеров дефектов обоих типов монотонно возрастает с увеличением энергии ПВА. При энергиях 4 и 5 КэВ, при которых наблюдается резкое снижение числа «выживающих» дефектов, образуются преимущественно одиночные дефекты и кластеры небольших размеров (до 10 вакансий и до 7 внедрений на кластер). В каскадах от ПВА энергий 15 и 20 КэВ увеличивается число дефектов, объединяющихся в кластеры довольно больших размеров. Так, при ЕПВА = 20 КэВ наблюдаются кластеры, содержащие до 118 внедрений, и кластеры, содержащие до 122 вакансий.
5 10 15
Энергия ПВА, КэВ
Рис. 5 Доля дефектов, образующих кластеры
"»Вшг
° Энергия ПВА, КэВ
Внедрения
Вакансии
Рис. 6 Число одиночных точечных дефектов и дефектов в кластерах, образующихся в каскадах смещений от ПВА различных энергии
Полученные нами результаты по кластерам дефектов качественно отличаются от результатов моделирования из работы [24], где при молекулярно-динамическом моделировании каскадов в железе, но с использованием другого многотельного потенциала и при температуре кристаллита Т=300 К, наблюдалось существенно большее образование кластеров внедрений, чем вакансий.
Заключение
В заключении кратко сформулируем основные результаты работы.
Путем компьютерного моделирования методом молекулярной динамики вычислены значения пороговых энергий смещения a-железа для различных кристаллографических направлений вылета ПВА. Величина средней пороговой энергии смещения составила 38,9 ± 1,5 эВ, что близко к рекомендуемому ASTM значению 40 эВ.
Исследована зависимость доли точечных дефектов, «выживающих» в каскаде смещений, от энергии ПВА (в диапазоне от 0,1 до 20 КэВ). Наименьшее значение доли «выживающих» дефектов наблюдается при ЕПВА = 10 КэВ и составляет =0,27. Обнаружено, что в диапазоне энергий ПВА от 3 до 5 КэВ число «выживающих» дефектов уменьшается с ростом энергии ПВА. Получены результаты по размерам и количеству кластеров вакансий и внедрений, образующихся в каскаде смещений, для различных энергий ПВА. Согласно этим результатам, количество вакансий, участвующих в процессе кластеризации, превосходит число образующих кластеры внедрений, и при ЕПВА > 5 КэВ доля дефектов обоих типов, образующих кластеры, возрастает с увеличением энергии ПВА.
Список литературы
1. Norgett, N. J. The proposed method of displacement doze rate calculation / N. J. Norgett, M. T. Robinson, I. M. Torrens // Nucl. Eng. And Design. - 1975. - V. 33. -Р. 50-56.
2. Wiedersich, H. Effects of the primary recoil spectrum on microstructural evolution / H. Wiedersich. - J. Nucl. Mater. -1991. - V. 70. - Р. 179-181.
3. Heinisch, H. L. Simulating the production of free defects in irradiated metals / H. L. Heinisch // Nuclear Instruments and Methods in Physics Research B 102. - 1995. -Р. 47-50.
4. Балашов, А. Н. Расчет методом молекулярной динамики первичных процессов радиационных повреждений, взаимодействия и свойств точечных дефектов в многокомпонентных материалах / А. Н. Балашов, М. Ю. Тихончев, Е. И. Шамари-на. - Препринт ТГТУ. - Тверь, 2000.
5. Daw, M. S. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals / M. S. Daw, M. I. Baskes // Phys. Rev. B 29. - 1984. -Р. 6443-6453.
6. Finnis, M. F. A simple empirical N-Body potential for transition metals / M. F. Fin-nis, J. E. Sinclair // Philos. Mag. A 50. - 1984. - Р. 45-55.
7. Rosato, V. Thermodynamical and structural-properties of FCC transition-metals using a simple tight-binding model / V. Rosato, M. Guellope, B. Legrand // Philos. Mag. A 59. - 1989. - № 2. - Р. 321-336.
8. Ackland, G. J. Development of an interatomic potential for phosphorus impurities in a-iron. / G. J. Ackland, M. I. Mendelev, D. J. Srolovitz // J. Phys. - 2004. - Condens. Matter 16. - Р. S2629-S2642.
9. Ackland, G. J. Computer simulation of Point Defect Properties in dilute Fe-Cu alloy using a many-body interatomic potential / G. J. Ackland, D. J. Bacon, A. F. Calder, T. Harry // Phil. Mag. A 75. - 1997. - Р. 713-732.
10. Pearson, W. B. A Handbook of Lattice Spacings and Structures of Metals and Alloys / W. B. Pearson. - New York : Pergamon, 1958.
11. Rayne, J. E. Elastic constants of iron from 4,2 to 300°K / J. E. Rayne, B. S. Chandrasekhar // Phys. Rev. - 1961. - V. 122. - Р. 1714-1716.
12. Schepper L. De, G. Knuyt, Stals L. M. [et al.] Materials Science Forum. -1987. - V. 15-18. - Р. 131-136.
13. Domain, C. Ab initio calculations of defects in Fe and dilute Fe-Cu alloys / C. Domain, C. S. Becquart // Phys. Rev. B 65. - 2002. - V. 024103. - Р. 024103-1-024103-14.
14. Kittel, C. Introduction to Solid State Physics / C. Kittel. - 6th ed. - New York : Wiley, 1986. - Р. 55.
15. Shewmon, P. G. Diffusion in Solids / P. G. Shewmon. - New York : McGraw-Hill, 1963.
16. Nordlund, K. Molecular dynamics simulations of threshold displacement energies in Fe. Nuclear Instruments and Methods / K. Nordlund, J. Wallenius, L. Malerba // Physics Research B. - 2006. - V. 246 (2). - Р. 322-332.
17. Vajda, P. Anisotropy of electron radiation damage in metal crystals / P. Vajda // Rev. Mod. Phys. - 1977. - V. 49. - № 3. - Р. 481-521.
18. Hohenstein, M. The anisotropy and temperature dependence of the threshold for radiation damage in gold - comparison with other FCC metals / M. Hohenstein, A. See-ger, W. Sigle // J. Nucl. Mater. - 1989. - V. 169. - Р. 33-46.
19. F. Maury, P. Anisotropy of the displacement energy in single crystals of molybdenum / Maury F., Vajda P., Biget M., Lucasson A., Lucasson P. // Radiat. Effects. 1975. -V. 25. - № 3. - Р. 175-185.
20. Maury, F. Anisotropy of defect creation in electron-irradiated iron crystals / F. Maury, M. Biget, P. Vajda, A. Lucasson, P. Lucasson // Phys. Rev. B 14. - 1976. -Issue 12. - Р. 5303-5313.
21. Lomer, J. N. Anisotropy of defect production in electron irradiated iron / J. N. Lomer, M. Pepper // Philos. Mag. - 1967. - V. 16. - Issue 144. - Р. 1119-1128.
22. Lucasson, P. G. Production and Recovery of Electron-Induced Radiation Damage in a Number of Metals / P. G. Lucasson, R. M. Walker // Phys. Rev. 127. - 1962. -Issue 2. - Р. 485-500.
23. ASTM E521, (E521-89) Practice for Neutron Radiation Damage Simulation by Charged-Particle Irradiation // Annual Book of ASTM Standards. - 1995. - V. 12.02.
24. Terentyev, D. A. Displacement cascades in Fe-Cr. A molecular dynamics study / D. A. Terentyev, L. Malerba, R. Chakarova [et al.] // J. Nucl. Mater. - 2006. -V. 349 (1). - Р. 119-132.
