Моделирование каскадов атомных смещений в бинарном сплаве Fe-9ат. %Cr методом молекулярной динамики Текст научной статьи по специальности «Физика»

УДК 621.039.531.001.57

МОДЕЛИРОВАНИЕ КАСКАДОВ АТОМНЫХ СМЕЩЕНИЙ В БИНАРНОМ СПЛАВЕ FE-9АТ. %CR МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ

© 2012 М.Ю. Тихончев, В.В. Светухин

Ульяновский государственный университет

Поступила в редакцию 20.11.2012

В работе представлены результаты компьютерного моделирования процессов первичной радиационной повреждаемости сплава Ре-9ат.%Сг методом молекулярной динамики. Моделирование проведено для температуры 300 К с использованием многотельных потенциалов межатомного взаимодействия. Рассмотрены каскады атомных смещений для энергий первично выбитого атома от 0.1 до 20 кэВ. Получены оценки количества дефектов, переживающих рекомбинацию в каскаде, а также результаты по количеству и размерам кластеров вакансий и собственных межузельных атомов, образующихся в таких каскадах. Установлено, что после прохождения каскада относительная хрома в СМА более чем вдвое превышает его концентрацию в исходном сплаве.

Ключевые слова: метод молекулярной динамики, каскад атомных смещений, точечный дефект, пара Френкеля.

1. ВВЕДЕНИЕ

Феррито-мартенситные стали с быстроспа-дающей наведенной активностью (RAFM) 7-10%-Cr-WVTa являются основными кандидатами для первой стенки и конструкций бланкета будущего реактора синтеза. Хотя RAFM-стали проявляют ряд преимуществ под облучением по сравнению с коммерческими сплавами, низкотемпературное упрочнение при нейтронном облучении сопровождаемое охрупчиванием, снижение ударной вязкости и пластичности остаются основным препятствием для их применения.

Данная работа посвящена моделированию каскадов атомных смещений в бинарном сплаве Fe-9ат.%Cr методом молекулярной динамики. Целью работы является определения ряда параметров первичной радиационной повреждаемости в указанном сплаве с учетом процессов рекомбинации и кластеризации точечных дефектов в каскадах смещений. Моделирование проведено для температуры 300 К.

Молекулярная динамика является одним из наиболее подходящих способов моделирования каскадов атомных смещений. К настоящему времени различными группами исследователе по всему миру проведено большое количество таких исследований применительно к различным материалам. Особенно много данных получено для

Тихончев Михаил Юрьевич, кандидат физико-математических наук, начальник лаборатории моделирования поведения неорганических материалов. E-mail: tikhonchev@sv.ulsu.ru

Светухин Вячеслав Викторович, доктор физико-математических наук, профессор, директор научно-исследовательского технологического университета. E-mail: slava@sv.uven.ru

чистого б-Бе с использованием различных потенциалов межатомного взаимодействия. Моделирование каскадов в сплавах также широко выполняется, хотя для многокомпонентных систем сохраняется проблема подготовки надежных потенциалов. В последние годы опубликовано достаточно большое количество работ, посвященных такому моделированию для системы Бе-Сг (см., например работы [1, 2, 3, 4, 5]). Малерба и др. [1] провели молекулярно динамическое моделирование каскадов смещений в сплаве Бе-10ат.%Сг для энергий первично-выбитого атома (ПВА) до 15 кэВ. Немного позднее Терентьев и др. [2] выполнили расчеты для энергий ПВА до 50 кэВ. Валлениус и др. [3] промоделировали каскады для системы Бе-5ат.%Сг и Бе-20ат.%Сг. В работах [1, 2, 3] использовались потенциалы межатомного взаимодействия на основе метода погруженного атома. Однако потенциал для Бе в работах [1, 2] показывает для собственных межузельных атомов (СМА) большую стабильность гантельной конфигурации <111>, чем <110>, что противоречит экспериментальным результатами и расчетам из первых принципов. Потенциалы в работе [3] дают почти вдвое более высокие энергии образования смешанных "гантелей" внедрений Бе-Сг в матрице а-Бе, чем показывают расчеты из первых принципов. Кроме перечисленного следует особо отметить, что потенциалы межатомного взаимодействия, использованные во всех упомянутых работах, предназначены для моделирования случайного сплава БеСг с заданной концентрацией Сг. Они не воспроизводят кривую энтальпии смешения Бе и Сг во всем диапазоне концентраций. Известно [6], что эта кривая меняет знак при концентрациях Сг ~5-

10%. В настоящей работе используется модифицированная версия потенциала А. Каро и др. [7], хорошо воспроизводящая энергию смешения случайного ферро-магнитного сплава БеСг для произвольной концентрации компонент.

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

Моделирование каскадов атомных смещений проводили методом молекулярной динамики. Мы использовали многотельный потенциал межатомного взаимодействия, предложенный А. Каро и др. [7]. В используемом формализме энергия системы из N атомов записывается как

Etot ^

1

TiTj ( x jWTTjir j)

(1)

(Ерщ (гУ ))+-X Ь

_ 1 * I 2 1 * I

где Eш - полная потенциальная энергия системы, Тг - тип г'-го атома: Бе или Сг, гг - собствен' ' Тг

ная электронная плотность атома Т как функция расстояния до его центра, Р - фунция внедрения, г. - расстояние между атомами г и _/,

УтгТг(г) - п'арный потенциал, УРеСг=(УРеРе+УСгСг)/2.

Н (х)=1, если Тг=Т)', иначе Н (х) - полином Редлиха-Кистера четвертой степени,

_ 1 Х _ 2

(

■ +

р ?л

pi р

(2)

j У

где pi - вклад в электронную плотность в узле i

/ „tot ч

( pi ) только от атомов железа.

Базовые потенциалы для Fe и Cr взяты авторами [7] из работ [8] и [3] соответственно. Параметры полинома h (x) были подобраны так, чтобы близко воспроизвести кривую энтальпии смешения случайного ферро-магнитного сплава Fe-Cr, рассчитанную методом ab initio из работы [6].

Мы произвели небольшую модификацию описанного потенциала. Во-первых, мы использовали более новую версию потенциала для Fe, взятую из работы [9]. Далее, необходимо аккуратно описать репульсивную часть потенциала для взаимодействий Fe-Cr при небольших (<1 Е) межатомных расстояниях. Это необходимо для моделирования каскадов атомных смещений [10,5]. Для этой цели мы заменили слагаемое hFeCr (x)VFeCr (r) в формуле (1) следующей функцией.

VPeCr (x, r ) _ h(x)VFC (Г ) + V(r ) ,

F'mod

■ PeCr

(r) _

Kcr(r),

r > r0

о

cr

+ dr2,0 < r < r

(3)

(4)

V(r) _ <

^л e _ , r . _

1 2 Ф(-), 0 < r < r1 4яе0 r a

a(r - r2 )3 e/r, r < r < r2

(5)

V(r) - экранированный кулоновский парный потенциал для r<r1. Длина экранирования a определена согласно работе [11], а функция F(x) -в соответствие подходу Циглера, Берсака и Лит-тмарка [12]. Коэффициенты c, d, а и // подбираются из условия непрерывности соответствующих функций и их первых производных, r=1 Е, rg и r2 - подгоночные параметры. Подгонка проведена под энергии связи смешанных FeCr гантелей внедрений с ориентацией <110> и <111>, рассчитанных Ольссоном и др. [13] методом ab initio. Полученные значения rg и r2 составили 2.22815 и 1.9664 Е соответственно.

3. МОДЕЛИРОВАНИЕ КАСКАДОВ СМЕЩЕНИЙ

Для моделирования каскадов атомных смещений и оценки числа "выживших" дефектов задавались ОЦК кристаллиты Fe-9ат.%Cr, содержащие до = 600000 атомов. При этом расположение атомов разных типов в сплаве носило случайный характер, т.е., каждому атому кристаллита с вероятностью p=0.91 приписывался тип Fe и с вероятностью (1-p) - Cr. При расчётах использовались "периодические" граничные условия. Поэтому, чтобы избежать цикличности возмущений вдоль плотноупакованного направления <111>, кристаллиты задавали в форме прямоугольного параллелепипеда, но не куба. Моделирование проводилось при начальной температуре кристаллита T = 600 K для восьми различных энергий ПВА (Епва): 0.1, 0.5, 1, 2, 5, 10, 15 и 20 кэВ. Начальную температуру обеспечивали путем задания начальных скоростей атомам кристаллита с последующим молекулярно динамическим моделированием NPT-ансамбля с баростатом и термостатом Берендсена для нулевого давления и температуры 300 К, соответственно, в течении 1 пс с шагом по времени 1 фс. Затем, фиксированному атому (ПВА) релаксированного кристаллита (это мог быть как атом Fe так и Ni) придавали импульс в высокоиндексном направлении <135>.

В ходе моделирования каскадов никакие алгоритмы сброса кинетической энергии для имитации "остывания" кристаллита не использовались. Рост температуры кристаллита составил от = 16 К для Епва=0.1 кэВ до = 125 К для ЕПВА=20 кэВ. Расчеты проводились с неравномерным шагом по времени, который выбирался так, чтобы он не превосходил 10-3 пс., и чтобы за один шаг по времени атом с максимальной кинетической энергией смещался не более чем на 0.02 Е. Моделируемое время развития каскада подбиралось так, чтобы обеспечить моделирование всего процесса образования и релаксации дефектов в каскаде вплоть до его затухания.

Fe

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

Доля "выживших" дефектов определялась по формуле:

Р (Е ПВА ) =

N ( Е ПВА )

(6)

/ ( Е ПВА )'

где М(Епва) - рассчитанное среднее число_"вы-живающих" дефектов, £(ЕПВА)=0.8-ЕПВА/(2 Еа ) -количество атомных смещений по МКТ-стандар-ту (без учета неупругих потерь энергии, которые в рассматриваемом интервале энергий ПВА являются незначительными).

Величину р(ЕПВА) в современной литературе часто называют коэффициентом каскадной эффективности. Для пороговой энергии смещения мы здесь использовали значение Еа =40 эВ. Полученные усредненные по направлениям импульса ПВА значения М(ЕПВА) и р(ЕПВА) представлены на рис. 1 и 2, соответственно, как для рассматриваемого сплава Ре-9ат.%Сг, так и для чистого а -Ре. Из полученных результатов видно, что везде на рассматриваемом интервале энергий наблюдается рост числа дефектов с увеличением Епва. Несмотря на то, что при ЕПВА. выше 15 кэВ наблюдается незначительное возрастание доли дефектов, переживших рекомбинацию в каскаде, с учетом погрешностей эту величину можно считать убывающей по ЕПВА.

Бэкон и др. [14] и Вудинг и др. [15] показали, что для металлов зависимость М(ЕПВА) хорошо аппроксимируется степенной функцией вида

К(Е) = А • Ев . (7)

Полученные нами результаты также хорошо описываются зависимостями такого вида (см. рис. 1), а именно

К(Е) = 5.45 • Е071, для Ре-9ат.%Сг, (8)

■ Ре-9ат.%Сг

о Риге Ре

Ре-9ат.%Сг аппроксимация Ре аппроксимация

1 10 Энергия ПВА, кэВ

Рис. 1. Рассчитанное число "выживающих" дефектов

■ Ре-9ат %Сг О Ре

-Ре-9ат %Сг аппроксимация

......Ре аппроксимация

0 5 10 15 20

Энергия ПВА, кэВ

Рис. 2. Доля дефектов, переживших рекомбинацию в каскаде смещений

К(Е) = 4.36 • Е068, для Ре, (9)

где Е - энергия ПВА в кэВ.

Учитывая (6) из формул (8-9) получаем соответствующие аппроксимации для зависимости Р(Епва):

К(Е) = 5.45 • Е~0'29, для Ре-9ат.%Сг, (10) р(Е) = 4.36 • Е~032, для Ре, (11)

На рис. 3 представлены доли хрома в СМА на момент завершения моделирования эволюции каскада. К этому моменту в кристаллите наблюдается лишь очень незначительная часть "гантелей" внедрений типа Сг-Сг, значительная часть межузельного хрома находится в смешанных "гантелях" Ре-Сг. Следует отметить, что средняя доля Сг в СМА слабо зависит от энергии ПВА и составляет = 22%. Таким образом, концентрация хрома в СМА более чем вдвое превышает его концентрацию в исходной матрице.

Известно, что существенный вклад в микроструктурную эволюцию материала под облучением вносит объединение производимых в нем точечных дефектов в кластеры. При моделировании каскадов смещений мы, наряду с оценкой числа "выживающих" дефектов, получили оцен-

10

0.1

0.1

1.2

0.8

5 0.6

0.4

0.2

0

40% 35% Й 30%

о

о

a io%

5%--

0% -I-,-,-,-,-1-,-,-,-,-1-,-,-,-,-1-,-,-,-,-1-,

0 5 10 15 20

Энергия ПВА, кэВ

Рис. 3. Доля хрома в СМА на момент затухания каскада

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

На рис. 4 представлены рассчитанные для сплава Бе-9ат.%Сг доли точечных дефектов, образовавших кластеры на момент завершения моделирования каскадов. Хорошо видно, что при энергиях ПВА 5 кэВ и выше количество вакансий, участвующих в процессе кластеризации, совпадает в пределах погрешности с числом СМА, входящих в кластеры.

4. ЗАКЛЮЧЕНИЕ

В работе методом молекулярной динамики проведено моделирование каскадов атомных смещений в бинарном сплаве Бе-9ат. %Сг для энергий ПВА от 100 эВ до 20 кэВ. При моделировании использовали многотельные межатомные потенциалы межатомного взаимодействия, взятые из работы [7]. Мы модифицировали эти потенциалы для корректного описания взаимодействий на небольших межатомных расстояниях, что важно для моделирования процессов радиационной повреждаемости.

Исследована зависимость доли дефектов, "выживающих" в каскаде смещений, от энергии ПВА (в диапазоне от 0.1 до 20 кэВ). Согласно полученным результатам зависимость числа "выживающих" дефектов хорошо аппроксимируется степенной функцией, что качественно согласуется результатами других исследователей. Значение каскадной эффективности при ЕПВА начиная с 10 кэВ составляет ~ 0.2. Эта оценка хорошо согласуется с результатами работы [5].

Энергия ПВА, кэВ

Рис. 4. Доля дефектов, образующих кластеры

В тоже время она лежит ниже соответствующей оценки ~ 0.3, полученной в работе [2], и превосходит оценку ~ 0.135 из работы [3]. В обоих случаях отличие может быть объяснено использованием различных потенциалов межатомного взаимодействия. В частности потенциал для Fe в работе [2] показывает для межузельных атомов большую стабильность гантельной конфигурации <111>, чем <110>, что противоречит экспериментальным результатами и ab initio расчетам для а -Fe. Кросс-потенциалы Fe-Cr в работе [3] дает заметно более высокие энергии образования смешанных "гантелей" внедрений Fe-Cr в матрице а-Fe.. Следует также отметить, что оценки как числа "выживающих" в каскаде пар Френкеля, так и каскадная эффективность для системы Fe-9ат.%Cr оказываются несколько выше, чем для чистого а-Fe.

Согласно полученным результатам доля хрома в СМА после прохождения каскада составляет = 22 %, т.е. более чем вдвое превышает долю хрома в исходном сплаве. Эти результаты качественно согласуются с результатами работы [2], но противоположны результатам, полученным в работах [4, 5]. Однако, несмотря на качественное согласие оценки доля Cr в СМА из работы [2] составляют 60 - 70 % для ЕПВА от 1 до 30 кэВ и около 50 % для Епва=50 кэВ. Эти отличия, по всей видимости, являются эффектом разных потенциалов. Для подтверждения или опровержения этих результатов необходимы дальнейшие, в том числе и экспериментальные, исследования.

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

Полученные результаты могут быть использованы для развития моделей радиационного повреждения сплавов Fe-Cr, в том числе, на основе многомасштабного подхода.

Работа выполнена при поддержке Минобрна-уки в рамках государственного задания на 20122014 гг, ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009 - 2013 годы и при частичной поддержке гранта РФФИ - проект 12-08-97076.

СПИСОК ЛИТЕРАТУРЫ

1. Molecular dynamics simulation of displacement cascades in Fe-Cr alloys / Malerba, L, Terentyev, D, Olsson, P., Chakarova R., Wallenius J. // J. Nucl. Mater. 329-333 (2004) p. 1156-1160.

2. Displacement cascades in Fe-Cr. A molecular dynamics study / Terentyev, D.A., Malerba, L, Chakarova, R, Nordlund, K, Olsson, P., Rieth, M, WalleniusJ. // J. Nucl. Mater., V. 349(1), 2006, p. 119-132

3. Modeling of chromium precipitation in Fe-Cr alloys / Wallenius J., Olsson P., Lagerstedt C., Sandberg N., Chakarova R., Pontikis V. // PHYSICAL REVIEW B 69, 094103 , 2004 094103-1-9

4. Lee HJ., Sadigh B., Shim J.-H., Wong K. Fundamental investigation of point defect interactions in Fe-Cr alloys. Presented at OECD SMINS Workshop Karlsruhe,

Germany, 6 June 2007.

5. MD simulation of atomic displacement cascades in Fe-10 at.%Cr binary alloy / M. Tikhonchev, V. Svetukhin, A. Kadochkin, E. Gaganidze //Journal of Nuclear Materials 395 (2009). P. 50-57.

6. Olsson P., Abrikosov I. A., Vitos L., and Wallenius J. // J. Nucl. Mater. 321, (2003) P. 84-90.

7. Caro A., Crowson D. A., and Caro M. // Phys. Rev. Lett. 95 (2005) 075702, 4 pp.

8. Mendelev M.I., Han S., Srolovitz D.J., Ackland G.J., Sun D.Y., andAsta M. // Phil. Mag. 83, 2003, p.3977-3994.

9. Development of an interatomic potential for phosphorus impurities in a -iron / G.J. Ackland, M.I.Mendelev, D.J.Srolovitz, S.W.Han, A.V.Barashev // J. Phys.: Condens. Matter 16 (2004) pp. S2629-S2642.

10. Giovanni Bonny and Lorenzo Malerba, OPEN REPORT of SCK^CEN-BLG-1022, 2005

11. Bohr N., Dansk Kgl.. Selsk Vid. // Mat. Fys. Medd., Vol. 18 (1948).

12. Ziegler J.F., Biersack J.P. and Littmark U. The stopping and range of ions in solids, Vol. 1 of 'The stopping and ranges of ions in matter', Pergamon Press, New York (1985).

13. Olsson P., Domain C., and WalleniusJ. // Phys. Rev. B 75, 2007, 014110, 12 pp.

14. Computer simulation of displacement cascade effects in metals / D.J. Bacon, A.F. Calder, F. Gao // Rad. Eff. Def. Sol. 141, 1997, p. 283 - 310.

15. A computer simulation study of displacement cascades in a -titanium / S.J. Wooding, D.J. Bacon, W.J. Phythian // Philos. Mag. A 72, 1995 p. 1261 - 1279.

SIMULATION OF ATOMIC DISPLACEMENT CASCADES IN BINARY FE^r %CR ALLOY BY MOLECULAR DYNAMICS METHOD

© 2012 M.Yu. Tikhonchev, V.V. Svetukhin

Ulyanovsk State University

This paper is devoted to the results of computer simulation of primary irradiation damage processes in Fe-9aT.%Cr by molecular dynamics method. Simulation is performed for temperature of 300 K. N-body interatomic potentials are used here. Cascade energies in rang of 0.1 - 20 keV are considered. Evaluations of amounts of defects survived in displacement cascades as well as numbers and sizes of vacancy and self-interstitial atoms (SIAs) are obtained. It is revealed that Cr fraction in SIAs is more than twice as higher in comparison with Cr fraction in initial alloy.

Key words: molecular dynamics method, atomic displacement cascade, point defect, Frenkel pair.

Mikhail Tikhonchev, Candidate of Physics and Mathematics, Head of Non-Organic Materials Behaviour Modelling Laboratory. E-mail: tikhonchev@sv.ulsu.ru Viacheslav Svetukhin, Doctor of Physics and Mathematics, Professor, Director of Research Institute of Technology. E-mail: slava@sv.uven.ru